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Introduction 


This  project  aims  to  improve  the  estimation  of  functional  properties  of 
breast  tissue  in  near  infrared  (NIR)  imaging  [1-3].  This  imaging 
technique  (also  known  as,  diffuse  optical  tomography  (DOT))  is  non- 
invasive  and  non-ionizing,  which  can  be  routinely  used  to  characterize 
the  breast  tissue.  In  this,  fibers  placed  on  the  boundary  of  breast  deliver 
NIR  light  (600  mn  -  950  nm)  and  collect  the  propagated  diffused  light 
(fiber-optic  setup  is  shown  in  Figure  1)  [2,3].  The  attenuation  and 
scattering  of  light  through  breast  tissue  volume  provide  an  estimation 
of  functional  properties  using  a  model-based  approach  [1].  This 
estimation  outcome  is  highly  dependent  on  the  model/algorithm,  more 
specifically  on  the  approach  to  match  the  experimental  data  with  the 
model  data  [1].  Even  though  the  light  propagates  in  three-dimensions 
(3D),  the  reconstruction  procedures  in  NIR  imaging  were  limited  to  two-dimensions  (2D)  due  to  computational 
complexity  and  limited  number  of  measurements.  In  this  context,  this  project  was  funded  to  explore  new 
reconstruction  methods  in  3D  and  decrease  the  computational  complexity  by  optimizing  these  procedures. 

Specific  aims  of  this  project,  in  brief 

1)  Reducing  the  computation  complexity  of  3D  imaging  by  investigating  different  data  collection  strategies 
and  optimizing  these  procedures. 

2)  Improving  the  quantitative  accuracy  of  optical  images  by  exploring  the  effect  of  penalty  terms  on  the 
reconstruction  techniques.  Incorporation  of  a  priori  information  from  other  modalities  (like  MRI,  CT) 
into  the  reconstruction  procedures  and  studying  its  effect. 

3)  Exploring  effective  ways  of  displaying  and  coregistering  3D  DOT  images. 

During  the  first  year  of  the  funding  period  of  this  project,  several  important  advances  towards  these  aims  have 
been  made.  Specifically,  optimization  of  critical  computational  aspects  in  NIR  imaging  was  completed.  An 
optimal  data-collection  strategy  especially  for  the  DOT  clinical  system  at  Dartmouth  for  the  current  estimation 
of  breast  tissue  optical  properties  was  also  found.  An  effective  way  for  usage  of  a-priori  structural  of 
information  from  MRI/CT  into  the  image  reconstruction  procedure  was  developed  and  proven  that  the 
quantitative  accuracy  of  DOT  images  can  be  improved  by  at  least  a  factor  of  two  with  this  additional 
information.  As  an  important  step  towards  realizing  the  final  outcome  of  this  project,  a  generalized  estimation 
procedure  was  developed  which  will  take  into  account  the  noise  characteristics  of  instruments  and  breast  tissue 
optical  properties  and  has  been  shown  robust  to  highly  noisy  data. 

Body 

Optimizing  the  critical  computational  aspects  of  near  infrared  tomographic  imaging 

The  image  resolution  and  contrast  in  Near-Infrared  (NIR)  tomographic  image  reconstruction  are  affected  by 
parameters  such  as  the  number  of  boundary  measurements,  the  mesh  resolution  in  the  forward  calculation  and 
the  reconstruction  basis.  The  magnitude  of  the  total  sensitivity  was  analyzed  to  find  the  spatial  variation  for  a 
given  problem,  and  the  field  response  of  the  domain  becomes  more  uniform  by  increasing  the  sensitivity  to 
deeper  regions,  while  suppressing  the  hypersensitivity  near  the  external  boundaries.  This  is  achieved  with  an 
increase  in  the  number  of  measurements. 

Using  singular-value  decomposition  (SVD)  and  example  reconstructed  images,  numbers  of  16  or  24  fibers  are 
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sufficient  for  imaging  the  2D  domain.  The  number  of  useful  measurements  actually  decreases  exponentially 
with  the  number  of  measurements  used,  and  the  number  of  useful  singular  values  increases  only  as  the 
logarithm  of  the  number  of  measurements.  For  this  2D  reconstruction  problem,  given  a  computational  limit  of 
10  sec  per  iteration,  leads  to  choice  of  forward  mesh  with  1785  nodes  and  reconstruction  (pixel)  basis  of  30x30 
elements. 


The  use  of  three  fundamentally  different 
data  collection  strategies  for  three- 
dimensional  (3D)  NIR  tomography  was 
compared.  Given  a  3D  NIR  imaging 
problem,  using  a  single  plane  of  data  can 
provide  useful  images  if  the  anomaly  to  be 
reconstructed  is  within  the  measurement 
plane.  However,  if  the  location  of  the 
anomaly  is  not  known,  3D  data  collection 
strategies  are  very  important.  The 
recovered  quantitative  accuracy  of  the 
anomaly  region  decreases  (approximately 
10%)  with  the  addition  of  out-of-plane  data 
relative  to  in-plane  data.  Usage  of  single¬ 
plane  of  data  gives  slightly  better 
quantitative  accuracy,  if  the  anomaly  lies  in 
the  data  acquisition  plane.  Further  the 
quantitative  accuracy  of  the  reconstructed 
anomaly  increased  approximately  from 
15%  to  89%  as  the  anomaly  moved  from 
the  centre  to  boundary,  respectively.  The 
data  supports  the  idea  that  the  use  of  in¬ 
plane  data  in  the  3D  data  collection 
strategies  may  be  sufficient  for  the  3D  NIR 
imaging. 

Complete  work  along  with  the  methods 
employed  and  detailed  discussion  of  the 
results  given  in  the  appendix  [4] 

(Yalavarthy  et  al,  Opt.  Express  14,  p.  6113-6127,  2006). 


3D  1-plane 


3D  31ayer:  inplane 


3D  3 layer:  out  of 
plane 


(a) 


(b) 


(c) 


0.008 


Figure  2:  The  reconstructed  absorption  coefficient  distribution  for 
the  cylindrical  object  with  a  spherical  absorption 
inhomogeneity  (diameter  of  15mm  and  contrast  2: 1  with  respect 
to  background)  located  at  x,  y  and  z  locations  (a)  (0,0,0), 

(b)  (30,0,0)  and  (c)  (30,0,10).  The  three  columns  of  images  show 
the  results  achieved  with  the  three  different  data  collection  schemes 


Usage  of  structural  a-priori  information 

NIR  tomography  combined  with  conventional  imaging  modalities  (MRI,  CT  and  Ultra  Sound)  has  been  a  very 
active  area  of  research.  These  hybrid  systems  show  superior  performance  in  terms  of  qualitative  (resolution) 
and  quantitative  accuracy  compared  to  stand-alone  systems.  But  still  there  is  lot  of  ambiguity  in  utilizing  the 
spatial  information  from  these  high  spatial  resolution  images  into  NIR  tomography  (coregistration).  This  work 
develops  a  simple  framework  to  incorporate  structural  a-priori  information.  Simple  weight  matrices  that  have 
Laplacian-type  or  Helmholtz-type  structures  that  are  derived  from  a-priori  information  have  been  developed.  It 
has  been  shown  that  utilization  of  structural  information  using  these  weight  matrices  will  not  bias  the 
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reconstruction  problem  towards  imperfect  structural  priors.  Usage  of  imperfect  a-priori  information  in  a 
parameter  reduction  (i.e.  hard-priors)  in  the  imaging  field  through  the  enforcement  of  spatially  explicit  regions 

gives  erroneous  results.  In  the  phantom 
experiments,  it  is  shown  that  the  Helmholtz  type 
of  regularization  matrix  gives  the  best  estimate  of 
the  scattering  parameter  and  the  Laplacian 
provides  best  estimate  for  the  absorption 
parameter.  Overall,  usage  of  structural-priors 
improve  the  reconstructed  image  quantitative 
accuracy  by  at  least  a  factor  of  two. 

Details  of  the  implementation  along  with  analysis 
of  results  given  in  the  appendix  [5]  (Yalavarthy  et 
al,  IEEE  Trans.  Med.  Imag.,  2006). 


Actual 


Laplacian  Helmholtz  Hard  Priors 


0  0  0  0 


0.006  0.008  0.01  0.012  0.014  0.016  0.018 


00 


0.02  0.04  0.06  0.08 


0.5  1 


0.6  0.8  1  1.2  1.4  1.6 

Figure  3:  Reconstruction  results  with  the  usage  of  imperfect  structural  priors  using 
different  reconstruction  techniques. 


Generalized  Least-Squares  (GLS)  minimization 


DOT  involves  recovery  of  the 
distribution  of  optical 
parameters  by  matching  the 
experimental  data  with 
modeled  data  (Levenberg- 
Marquardt  (LM) 

Minimization).  A  variation  of 
this  approach  by  adding  the 
parameter  field  to  the 
minimization  function  is  done 
through  Tikhonov 

regularization,  where  the 
regularization  parameter  is 
chosen  to  overcome  the  ill- 
conditioning  of  the  problem.  In 
this  work,  a  generalized 
framework  for  DOT  was 
developed  including  variance 
of  the  data  and  parameters  as  weight  matrices.  These  weight  matrices  can  also  include  the  structural  information 
obtained  by  MRI,  Ultrasound  or  X-ray  imaging.  These  weight  matrices,  include  the  system  noise  characteristics 
and  expected  size  of  optical  parameters  and  constraints  for  the  imaging  problem  and  make  the  inversion  routine 
more  robust  to  noise.  This  also  makes  the  imaging  problem  more  stable.  It  is  also  important  to  note  that 
Tikhonov  regularization  becomes  a  special  case  of  the  Generalized  Least-Squares  (GLS)  formulation.  This  GLS 
estimation  of  optical  properties  has  been  shown  to  be  very  robust  to  noise  and  proven  to  be  stable  over 
iterations. 


Actual  LM  Tikhonov  GLS  AC 


0.009  0.01  0.011  0.012  0.013  0.014  0.015  0.016  0.017  0.016  0.019 


1  1.2  1.4  1.6  2 

Figure  4:  Reconstruction  results  using  different  minimization  techniques  with  3% 
noise  in  the  data. 
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Complete  formulation  along  with  results  is  given  in  the  appendix  [6]  (Yalavarthy  et  al,  Med.  Physics,  2007). 

Key  research  accomplishments 

•  Optimization  of  computational  aspects  of  DOT  image  reconstruction  (especially,  data-collection 
strategy) 

•  Effective  way  of  incorporating  structural  priors  into  NIR-DOT  image  reconstruction  procedure 

•  Generalized  least  squares  minimization  formulation  and  extensive  testing  of  the  algorithm  in  simulations 

Reportable  outcomes 

First  year  of  this  training  program  has  lead  to  two  peer-reviewed  journal  publications  and  a  number  of 
conference  presentations  (and  proceedings  papers  as  well).  In  detail: 

1.  Poster  presentation  at  Optical  Society  of  America  Biomedical  Optics  Topical  Meeting,  Florida,  March 
19-22,2006. 

2.  Poster  presentation  at  Network  for  Translational  Research  Optical  Imaging  Network  (NTROI)  Retreat, 
Hyatt  Regency  Newport  Beach,  CA,  June  22-24,  2006. 

3.  Oral  presentation  the  International  Society  for  Optical  Engineering  (SPIE)  BiOS-2007  in  Photonics 
West-2007,  San  Jose,  California,  20-25  January  2007. 

Conclusions 


This  project  is  part  of  continuing  effort  to  develop  methods/algorithms  for  three-dimensional  alternative  breast 
imaging  modalities  at  Dartmouth.  Some  important  miles  stones  in  the  project  include  completing  the  work  on 
optimizing  the  NIR  data-collection  strategies  in  3D  (completing  Aim-1).  A  framework  to  incorporate  the 
spatial-priors  in  to  the  NIR  image  reconstruction  procedure  was  developed  and  also  proven  to  be  effective  even 
in  case  of  imperfect  spatial  priors,  which  is  part  of  Aim-2.  Moreover,  a  new  algorithm  that  takes  into  account 
noise  characteristics  of  the  instruments  was  developed  and  tested  extensively  in  the  simulation  studies. 
Preliminary  3D  reconstruction  results  using  this  new  algorithm  show  improved  quantitative  accuracy  compared 
to  the  traditional  image  reconstruction  techniques.  Finally  steps  are  taken  towards  parallelizing  the  codes 
developed  here  to  reduce  the  run  time  and  memory  requirements. 
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Appendices 

Three  journal  papers  (manuscripts)  are  attached  which  were  a  direct  result  of  this  project. 
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Abstract:  The  image  resolution  and  contrast  in  Near-Infrared  (NIR) 
tomographic  image  reconstruction  are  affected  by  parameters  such  as  the 
number  of  boundary  measurements,  the  mesh  resolution  in  the  forward 
calculation  and  the  reconstruction  basis.  Increasing  the  number  of 
measurements  tends  to  make  the  sensitivity  of  the  domain  more  uniform 
reducing  the  hypersensitivity  at  the  boundary.  Using  singular-value 
decomposition  (SVD)  and  reconstructed  images,  it  is  shown  that  the 
numbers  of  16  or  24  fibers  are  sufficient  for  imaging  the  2D  circular  domain 
for  the  case  of  1%  noise  in  the  data.  The  number  of  useful  singular  values 
increases  as  the  logarithm  of  the  number  of  measurements.  For  this  2D 
reconstruction  problem,  given  a  computational  limit  of  10  sec  per  iteration, 
leads  to  choice  of  forward  mesh  with  1785  nodes  and  reconstruction  basis  of 
30x30  elements.  In  a  three-dimensional  (3D)  NIR  imaging  problem,  using  a 
single  plane  of  data  can  provide  useful  images  if  the  anomaly  to  be 
reconstructed  is  within  the  measurement  plane.  However,  if  the  location  of 
the  anomaly  is  not  known,  3D  data  collection  strategies  are  very  important. 
Further  the  quantitative  accuracy  of  the  reconstructed  anomaly  increased 
approximately  from  15%  to  89%  as  the  anomaly  is  moved  from  the  centre 
to  boundary,  respectively.  The  data  supports  the  exclusion  of  out  of  plane 
measurements  may  be  valid  for  3D  NIR  imaging. 
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1.  Introduction 

In  the  recent  years,  there  has  been  a  heightened  interest  in  near- infra-red  (NIR)  optical 
tomography,  for  applications  such  as  diagnostic  breast  cancer  imaging  [1-3]  and  for  brain 
function  assay  [1,  4,  5].  In  NIR  tomography,  the  aim  is  to  reconstruct  interior  optical 
properties  of  the  tissue  under  investigation  from  a  finite,  yet  incomplete  set  of  transmission 
measurements  taken  at  the  tissue  external  boundaries.  The  reconstructed  optical  properties  can 
give  clinically  useful  information  regarding  tissue  physiology  and  state,  such  as  chromophore 
concentration  and  oxygen  saturation.  Typically,  the  optical  source  light  used  for  excitation  in 
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NIR  studies  is  delivered  through  optical  fibers  and  the  transmitted  light  is  also  collected 
through  the  same  or  additional  fibers  which  are  in  contact  with  the  external  surface  of  the 
tissue.  Using  these  measurements,  distributions  of  wavelength  dependent  absorption  and/or 
scattering  coefficients  of  the  tissue  are  reconstructed  using  a  model-based  iterative  algorithm. 
NIR  studies  have  the  advantage  of  being  non-invasive,  non-hazardous  and  can  therefore  be 
applied  repeatedly  to  investigate  functional  changes  in  tissue  over  a  prolonged  time. 

The  dominance  of  light  scattering  in  tissue  at  NIR  wavelengths  makes  optical  tomography 
inherently  more  difficult  in  the  sense  that  light  becomes  diffuse  within  millimeters  of  travel, 
reducing  the  resolution  of  the  reconstructed  images.  The  image  reconstruction  procedure  (i.e. 
the  inverse  problem)  is  non-linear,  ill-posed  and  ill-conditioned  [6]  and  to  improve  image 
reconstruction,  the  number  of  measurements  are  generally  increased,  to  increase  the  amount 
of  independent  information.  However  due  to  experimental  set-up  constraints,  such  as  the  light 
collection  strategy,  source  and  detector  fiber  size  and  the  imaging  domain  geometry,  the  total 
number  of  boundary  measurements  that  can  be  taken  from  is  often  quite  limited.  In  addition, 
there  are  constraints  on  the  data  acquisition  and  computation  time  that  need  to  be  considered 
for  the  specific  application  in  which  NIR  light  is  used. 

There  have  been  some  limited  studies  [7-11]  on  optimization  of  the  fiber  positions  and 
measurements  to  get  the  best  possible  image  resolution  and  contrast  in  NIR  tomography. 
More  specifically.  Culver  et.  al  [11]  have  showed  that  singular  value  decomposition  (SVD) 
analysis  of  the  weight  matrix  (also  known  as  the  Jacobian  or  sensitivity  matrix)  can  be  used  to 
optimize  detector  placement  in  the  reflectance  and  direct  transmittance  geometries  of  a 
homogeneous  slab  medium,  and  indicated  that  this  could  be  extended  to  arbitrary  geometries 
with  heterogeneous  tissue  volumes.  However,  there  remain  many  unknowns  regarding  the 
appropriate  number  of  measurements  required  to  get  a  sufficiently  good  image  given  the 
practical  constraints  of  measurement  number  and  image  recovery  algorithm,  which  is  the 
subject  of  this  paper.  Furthermore,  few  studies  have  specifically  investigated  the  effect  of 
mesh  resolution  in  both  the  forward  and  inverse  calculations  and  very  little  is  known  about  the 
quantitative  increase  in  accuracy  which  is  a  direct  result  of  mesh  resolution  and  appropriate 
reconstruction  bases.  This  work  is  an  attempt  to  answer  questions  regarding  the  limited 
increase  in  number  of  measurements,  more  specifically  benefits  from  the  increased  amount  of 
information  as  well  as  investigating  aspects  that  will  have  effects  on  image  reconstruction 
procedure  and  resolution  as  well  as  the  contrast  of  the  reconstructed  image. 

In  the  present  work,  both  a  two  dimensional  (2D)  circular  domain  and  a  three  dimensional 
(3D)  cylindrical  geometry  are  investigated  since  most  investigations  to  date  have  used  either 
of  these  geometries  for  system  and  algorithm  evaluation.  Initially  the  effect  of  mesh  resolution 
is  investigated  in  the  forward  problem  by  comparing  the  Jacobian  cross-section  for  various 
resolution  2D  meshes  to  show  improvements  in  numerical  accuracy.  Next  the  effect  of 
increasing  the  number  of  measurements  upon  the  resulting  reconstructed  image  using 
singular-value  analysis  is  investigated.  Results  regarding  the  optimized  reconstruction  basis 
are  presented  for  the  given  2D  model,  and  the  impact  in  the  Root  Mean  Square  (RMS )  error  of 
increased  spatial  sensitivity  is  presented  as  a  function  of  increasing  number  of  measurements. 
Finally  a  case-to-case  analysis  is  shown  by  increasing  the  number  of  measurements  in  image 
reconstruction  procedure  and  comparing  the  underlying  image  errors  within  the  reconstructed 
images. 

Since  3D  problems  have  more  degrees  of  freedom  (unknown  parameters),  they  are  highly 
ill-determined  as  compared  to  the  2D  problem.  But  NIR  optical  tomography  utilizes  the  data 
from  the  3D  tissue  volumes  and  therefore  should  be  treated  as  a  3D  imaging  problem.  Since 
light  propagation  in  tissue  is  physically  spread  in  all  directions,  3D  models  are  known  to  be  an 
accurate  prediction  of  the  light  fluence,  whereas  2D  models  are  simple  yet  inaccurate  at 
predicting  the  interior  fluence  distributions  [4,  12-17].  In  order  to  further  advance  NIR  optical 
tomography  into  a  suitable  and  accurate  clinical  imaging  modality,  it  is  important  to  develop 
fully  3D  imaging  tools,  yet,  the  major  challenge  in  this  task  is  to  determine  how  to  acquire 
large  data  sets  which  overcome  the  inherent  limitation  of  the  3D  problem  being  ill-determined 
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[18].  That  is,  to  improve  image  reconstruction  quality  in  3D,  the  number  of  measurements  can 
be  increased  as  mentioned  in  2D  case,  even  here  these  measurements  are  quite  limited. 

For  the  chosen  3D  cylindrical  geometry,  for  example,  acquiring  experimental  data  from 
three  different  planes  of  fiber  setup  improves  the  reconstructed  image  of  the  entire  domain  as 
compared  to  one  single  plane  of  data,  as  there  are  greater  numbers  of  measurements  providing 
a  larger  set  of  sampling  of  the  entire  volume  of  interest.  There  are  many  strategies  to  increase 
measurement  number  and  it  is  not  clear  which  present  the  best  improvement  in  the  final 
image.  Specifically,  this  work  examines  effects  of  different  measurement  strategies  for  3D 
NIR  tomography  by  presenting  and  quantifying  the  underlying  effects  of  using  a  single  plane 
of  tomographic  data  as  compared  to  three  planes  of  tomographic  data.  Within  the  latter  case, 
this  work  also  presents,  quantifies  and  discusses  the  benefits,  limits  and  losses  due  to  the 
measurement  of  in-plane  data  as  compared  to  out-of  plane  data  and  will  compare  and  contrast 
these  data  collection  geometries  from  the  prospective  of  gain  and  loss  in  the  reconstructed 
image  quality  and  respective  computation  time. 

2.  Methods 


Conventional  numerical  methods  for  the  forward  calculations  in  NIR  imaging  use  the  Finite 
Element  Method  (FEM),  which  is  considered  as  a  flexible  and  accurate  approach  to  modeling 
heterogeneous  domains  with  arbitrary  boundaries.  Light  transport  in  scattering  tissue  can  be 
accurately  described  by  the  Diffusion  Approximation  (DA)  to  the  Radiative  Transfer  Equation 
(RTE)  [19]: 


-  V.r(r)Vi>(r,ffl)  + 


A„(D- 


<S>(r,0})  =  q0(r,0)) 


(1) 


where  0(r,<w)is  the  photon  density  at  position  r  and  modulation  frequency  CO  (100  MFlz  in 

this  work),  and  K  =  l/[3(|Ua  +  |4S )],  the  diffusion  coefficient,  where  |4a  and  (Xs  are  the 
probabilities  per  unit  length  of  absorption  and  transport  scattering,  respectively,  and 

<70(r,  CO)  is  an  isotropic  source  term.  The  Robin  (Type  III)  boundary  condition  is  used  which 

best  describes  the  light  interaction  from  a  scattering  medium  to  the  external  air  boundary  [20]. 
The  calculated  boundary  data  values  with  a  frequency  domain  system  are  the  amplitude  and 
phase  of  the  signal,  from  which  the  diffusion  and  absorption  coefficients  can  be 
simultaneously  reconstructed. 

For  the  inverse  problem,  a  small  change  in  boundary  data  is  related  to  a  small  change  in 
optical  properties  through  the  Jacobian  matrix  of  values.  The  Jacobian  matrix  for 
reconstructing  both  the  unknowns  using  two  different  data-types  is  calculated  using  the 
Adjoint-method  [21],  and  has  dimensions  of  (2xSxD)  by  (2xN),  where  S  and  D  are  the 
number  of  sources  and  detectors  corresponding  to  each  source  respectively.  N  represents  the 
number  of  nodes  in  the  mesh  used  in  the  forward  calculation.  Here  the  Jacobian  maps  the 
changes  in  log  amplitude  and  phase  (2xSxD)  to  both  absorption  and  diffusion  changes  at  each 
node  of  the  FEM  model  (2xN).  The  Jacobian  which  maps  the  change  in  detected  signal  to 
image  space  has  four  parts: 


j  = 


die  dfia 
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In  all  our  analysis,  only  the  72  section  is  considered  (dimension  of  (SxD)  by  N),  which  maps  a 
small  change  in  the  absorption  coefficient  to  a  small  change  in  measured  log  intensity  of  the 
signal.  Since  all  kernels  of  the  complete  Jacobian  show  similar  results,  the  discussion  is 
limited  to  the  results  of  72,  and  shall  henceforth  be  referred  to  as  J. 

In  the  reconstruction  procedure  presented,  a  modified  Levenberg-Marquadt  algorithm  is 
used  for  calculating  the  estimates  of  jaa,  which  is  an  iterative  procedure  [10]  solving: 
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[And  =  [JTJ  +  XI]'1.  JTb 


(3) 


Here  [  A|j,a  |  is  an  update  vector  for  the  absorption  coefficient,  I  is  the  identity  matrix  and  X  is  a 
regularization  parameter.  Also,  b  =  [y  -  jFYjla)L  where  y  is  the  measured  (or  simulated) 
heterogeneous  boundary  data  and  JXl-fi )  is  the  forward  data  for  the  current  estimate  of  |4a.  In 
all  of  the  presented  work  using  simulated  data,  1%  noise  was  added  to  the  amplitude,  which  is 
a  typical  noise  observed  in  experimental  data  [2], 

For  the  2D  analysis  a  circular  model  with  a  diameter  of  86  mm  centered  at  (0,  0)  and  with 
homogeneous  optical  properties  of  |4a  =  0.01  mm'1  and  ps  =  1.0  mm"1  is  considered.  The  light 
collection/delivery  fibers  are  arranged  in  a  circular  equally  spaced  fashion,  where  one  fiber  is 
used  as  the  source  while  all  other  fibers  are  used  as  detectors,  to  give  ‘P’  number  of 
measurements  (where  P=  M(M-l),  where  M  is  number  of  fibers).  The  source  is  a  Gaussian 
source  of  Full  Width  Half  Maximum  (FWHM)  of  3mm,  and  it  is  placed  one  transport 
scattering  length  within  the  external  boundary. 


Fig.  1 .  Schematic  diagram  of  data  collection  geometry  used  for  the  3D  cylindrical  model. 

For  the  3D  analysis,  a  cylindrical  medium  with  a  diameter  of  86  mm  having  height  of  100 
mm  centered  at  (0,  0,  0),  with  homogeneous  optical  properties  of  |4a  =  0.01  mm"1  and  (i/  =  1 
mm"1  is  used  (Fig.  1).  The  light  collection/delivery  fibers  are  arranged  in  a  circular  and 
equally  spaced  fashion  and  are  in  either  a  single  plane  of  16  fibers  or  3  planes  of  16  fibers  per 
plane,  totaling  48  fibers.  Specifically  three  different  strategies  for  data  collection  are 
considered: 

(a) .  Single  layer  data:  The  16  fibers  are  arranged  in  a  circular  and  equally  spaced  fashion  in  a 
single  Layer-I  (Fig.  1 ),  where  one  fiber  is  used  at  a  time  as  the  source  while  all  other  fibers  are 
used  as  detectors,  to  give  240  ( 16x15)  amplitude  measurements. 

(b) .  Three  layers  of  in-plane  data:  The  48  fibers  are  arranged  in  a  circular  equally  spaced 
fashion  in  all  three  layers  (Layers-I,  II  &  III  in  Fig.  1),  giving  16  fibers  per  plane,  where  one 
fiber  is  used  at  a  time  as  the  source  while  only  those  fibers  in  the  same  “source  fiber  layer”  are 
used  as  detectors,  to  give  720  (3x16x15)  amplitude  measurements. 

(c) .  Three  layers  of  out-of-plane:  Same  as  above,  except  when  one  fiber  is  used  at  a  time  as 
the  source,  all  other  fibers  in  all  three  planes  are  used  as  detectors.  This  leads  to  2256  (48x47) 
amplitude  measurements. 

For  the  image  reconstruction  process,  an  iterative  update  to  the  Jacobian  matrix  was 
computed,  after  each  successive  image  estimation.  At  each  iteration,  the  objective  function 
was  evaluated  to  estimate  the  projection  error.  The  reconstruction  procedure  was  then  stopped 
when  the  projection  error  decreased  by  less  than  3%. 
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Fig.  2.  The  sensitivity  (Jacobian)  contour  plot  of  log  amplitude  and  |ia  for  a  source  (S)  and 
detector  (D),  which  are  diagonally  opposite  to  each  other  as  shown,  calculated  on  a  circular 
mesh  of  9664  nodes. 

2.1.  2D  Mesh  Resolution 

In  FEM  the  domain  is  divided  into  finite  discretized  sub-domains  wherein  the  numerical 
accuracy  and  stability  depends  highly  on  this  discretization  (mesh  resolution).  Since  the 
Jacobian  represents  the  sensitivity  of  the  detected  signal  to  a  small  change  in  optical 
properties,  the  numerical  accuracy  of  this  value  is  crucial  component  of  the  image 
reconstruction  problem,  to  study  the  effect  of  mesh  resolution  in  2D  case,  we  choose  different 
resolution  meshes  (with  number  nodes  ranging  from  150  to  4617  nodes)  along  with  a  high- 
resolution  mesh  of  9664  nodes  for  calculation  of  Jacobian.  The  Jacobian  with  a  diagonally 
opposite  source  and  detector  is  used,  as  shown  in  Fig.  2,  from  which  the  RMS  error  is 
calculated  for  each  mesh  with  respect  to  the  high-resolution  mesh.  The  RMS  error  is 
calculated  by  interpolating  the  Jacobian  of  each  mesh  unto  a  uniformly  distributed  grid, 
allowing  direct  comparison  of  each  result.  Since  the  Jacobian  represents  the  sensitivity  of  the 
detected  signal  to  a  small  change  in  optical  properties,  the  numerical  accuracy  of  this  value  is 
a  crucial  component  of  the  image  reconstruction  problem.  Here  the  highest  resolution  mesh 
provides  the  most  accurate  and  numerically  stable  solution,  therefore  the  calculated  RMS 
error  indicates  the  numerical  accuracy  of  each  lower  resolution  mesh.  The  computation  time 
taken  for  calculation  of  Jacobian  and  forward  data  is  also  noted  as  a  function  of  mesh 
resolution.  All  the  computations  were  carried  out  on  Pentium-IV  2.5  GHz  processor  with  2 
GB  of  RAM. 

2.2.  Singular-Value  (SV)  analysis 

Singular- Value  (SV)  analysis  for  the  Jacobian  matrix  is  explained  in  detail  elsewhere  [10]. 
Using  SV-analysis,  the  Jacobian  is  decomposed  into: 

J  =  USVT  (4) 

where,  U  &  V  are  orthonormal  matrices  containing  the  eigenvectors  of  J  and  S  is  a  diagonal 
matrix  containing  the  singular  values  of  J.  Vectors  of  U  and  V  correspond  to  the  modes  in  the 
detection  space  and  image  space,  respectively,  while  the  magnitude  of  the  singular  values  in  S 
represents  the  importance  of  the  corresponding  eigenvectors  in  U  and  V.  More  nonzero 
singular  values  indicating  more  modes  are  effective  in  between  the  two  spaces,  which  bring 
more  detail  and  improve  the  resolution  in  the  space.  There  are  normally  P  nonzero  singular 
values  in  the  diagonal  matrix  and  these  values  are  sorted  in  decreasing  order.  Typically  only 
those  singular  values  above  the  noise  level  (in  this  study,  1  %  noise  in  amplitude)  are  used,  as 
they  contain  the  only  useful  information  in  the  matrix.  Thus,  it  is  possible  to  detennine 
whether  increasing  the  number  of  measurements  gives  rise  to  an  increase  in  the  number  of 
useful  singular  values,  which  indicates  improvement  in  the  recovered  images. 
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In  2D,  this  analysis  was  applied  to  two  separate  cases:  (1)  a  homogeneous  case  with  optical 
properties  as  given  before,  and  (2)  a  heterogeneous  case  which  mimics  breast  optical 
properties  [22],  with  properties  of  fibro-glandular  layer  being  pa  =  0.003  mm'1  and  ji/  =  0.95 
mm'1  and  having  diameter  of  66  mm  and  fatty  layer  surrounding  it  having  jua  =  0.006  mm'1 
and  jj./  =  1.1  mm'1  with  a  thickness  of  20mm.  The  number  of  useful  singular  values  above  the 
noise  level  were  calculated  as  the  number  of  measurements  was  increased.  The  mesh  that  was 
found  to  have  an  optimum  resolution  from  the  previous  analysis  of  the  Jacobian  (Sec.  2. 1 )  was 
used  for  these  analysis.  For  both  these  cases,  the  percentage  of  useful  measurements  with 
respect  to  total  number  of  measurements  was  calculated  as: 


Useful  measurements  (in  %)  = 


Useful  number  of  singular  values 
Total  number  of  singular  values 


(5) 


Additionally,  the  effect  of  mesh  resolution  was  studied  for  its  impact  on  the  number  of 
independent  boundary  data  points  with  an  increase  in  number  of  measurements  by  calculating 
the  rank  of  the  Jacobian,  which  is  defined  as  the  maximum  number  of  linearly  independent 
rows/columns  of  a  given  matrix.  As  each  row  of  the  Jacobian  indicates  each  measurement,  the 
rank  of  the  Jacobian  indicates  the  total  number  of  independent  measurements. 

Image  reconstruction  consists  of  two  separate,  yet  equally  important  parts;  the  forward 
model  and  the  inverse  model.  For  the  forward  model,  the  mesh  used  in  FEM  needs  to  be  such 
that  to  ensure  numerical  accuracy,  as  already  discussed.  For  the  inverse  problem,  however,  the 
goal  is  to  reduce  the  number  of  unknowns  for  the  iterative  update  by  the  use  of  a 
reconstruction  basis  [23].  Therefore  it  is  important  to  investigate  the  effects  of  various 
reconstruction  basis  degrees  of  freedom  on  the  reconstruction.  Various  reconstruction  basis 
can  be  used,  such  as  second  mesh  basis  [24],  pixel  basis  [23]  or  adaptive  [25,  26]  .  With  this 
goal,  a  reconstruction  basis  was  optimized  for  the  given  2D  problem  by  looking  at  the  number 
of  useful  singular  values  for  various  pixel  (reconstruction)  basis.  A  linear  pixel  basis  of  having 
100  (10  by  10)  elements  to  1600  (40  by  40)  elements  was  used  and  the  Jacobian  was  mapped 
to  this  basis  for  the  analysis. 


Table  1.  The  RMS  error  (with  respect  to  the  fine  mesh  of  9664  nodes)  in  the  Jacobian  cross-section  from  center  to 
boundary,  (indicated  by  dashed  line  in  Fig.  2)  at  y  =  0  mm.  This  is  tabulated  as  a  function  of  mesh  resolution,  or 
number  of  nodes  in  the  mesh.  Last  two  rows  show  the  computation  time  taken  for  calculation  the  Jacobian  and 
Forward  data  for  16  source- detector  pairs  (240  measurements).  For  the  fine  mesh  of  9664  nodes  the  computation  time 
for  Jacobian  and  Forward  data  is  98.1  sec  and  28  sec  respectively. 


Nodes 

150 

425 

1360 

1785 

2683 

3047 

3569 

4617 

RMS  error 

60.56 

27.84 

5.06 

4.84 

2.57 

2.15 

1.85 

1.07 

Jacobian 
Computation 
Time  (in  Sec.) 

1.1 

2.5 

7.8 

10.1 

15.2 

17.8 

20.8 

38.1 

Forward  data 
Computation 
Time  (in  Sec.) 

0.1 

0.3 

0.9 

1.2 

1.9 

2.2 

2.6 

9.8 

2.3.  Reconstruction  examples 


In  order  to  understand  the  effect  of  increasing  the  number  of  measurements  on  total  sensitivity 
for  a  given  2D  model  the  magnitude  of  the  Jacobian  was  examined  as  a  function  of  number  of 
measurements.  To  achieve  this,  the  horizontal  cross-section  of  the  whole  Jacobian  was 
plotted,  which  was  summed  up  over  all  measurements,  from  center  to  boundary,  and 
examined  as  the  number  of  measurements  increased.  Since  the  Jacobian  provides  relative 
sensitivity,  a  cross-section  plot  was  normalized  in  each  case  with  respect  to  its  magnitude  at 
the  center  of  the  model  and  calculated  as  a  function  of  number  of  measurements  (56  to  4032). 
For  the  3D  model,  the  cross-section  of  the  total  Jacobian  was  normalized  with  respect  to  its 
magnitude  at  the  center  of  the  model  (as  in  the  2D  case),  for  each  case  of  the  three  different 
data  collection  strategies.  Finally,  for  the  2D  model,  only  the  absorption  coefficient  was 
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reconstructed  with  an  increasing  number  of  measurements  of  an  object  with  absorption 
inhomogeneity  at  various  positions  of  domain  using  log  of  amplitude  data.  A  circular 
absorption  anomaly  of  diameter  of  10  mm  was  used  having  a  contrast  of  2:1  compared  to  its 
background.  We  used  the  optimal  forward  mesh  along  with  optimal  reconstruction  basis  for 
the  reconstruction  procedure.  A  total  of  2  positions  of  absorption  inhomogeneity  were 
considered  with  it  center  at  (x,y)  of  (0,  0),  and  (30,  0)  for  various  number  of  measurements 
starting  from  56  to  4032. 


Fig.  3.  Singular  value  analysis  of  homogeneous  and  heterogeneous  2D  circular  models,  (a). 
Plot  of  the  useful  singular  values  versus  number  of  measurements,  (b).  Plot  of  percentage  of 
useful  measurements  versus  the  total  number  of  measurements,  (c).  Plot  of  the  Rank  versus 
number  of  measurements  is  shown  for  a  range  of  mesh  nodes,  (d).  Plot  of  the  number  of  useful 
singular  values  versus  number  of  measurements  is  shown,  for  various  reconstruction  bases. 


For  the  3D  case,  a  spherical  absorption  anomaly  of  diameter  of  15  mm  was  assumed 
having  a  contrast  of  2:1  compared  to  its  background.  A  total  of  3  positions  of  absorption 
inhomogeneity  were  considered  with  its  center  at  x,  y  and  z  of  (0,0,0),  (30,0,0)  and  (30,0,10). 
The  anomalies  were  reconstructed  using  the  noise  added  data  (1%  in  amplitude)  simulated 
from  the  three  different  fiber  location  strategies.  Full  Width  at  Half -Maximum  (FWHM)  was 
measured  for  each  of  the  peaks  in  the  X-Y  and  Z-Y  planes  as  well  as  the  total  computation 
time  for  reconstruction  process. 


Table  2.  The  number  of  useful  measurements  above  the  1%  expected  noise  level,  is  shown  for  the  2D  circular  and  3D 
cylindrical  models,  having  16  source  and  detector  fibers  with  one  or  three  planes  of  data  collection.  The  two  upper 
rows  have  only  1  plane  of  collection,  whereas  the  2nd  last  row  has  3  planes  of  collection  but  not  between  the  planes, 
and  the  last  row  has  3  planes  of  data  collection  with  complete  out  of  plane  measurements. 


Number  of 
Unknowns 

Number  of 
Measureme 
nts 

Number  of 
Useful  Singular 
values 

Useful 

measurements 

(%) 

Magnitude  of  largest 
singular  value 

2D 

1785 

240 

91 

37.92 

796.4 

3D  1  layer 

20163 

240 

107 

44.58 

117.1 

3D  3layer  in¬ 
plane 

20163 

720 

269 

37.36 

164 

3D  3layer  out- 
of-plane 

20163 

2256 

328 

14.54 

304.6 

3  Results 

Figure  2  shows  a  sensitivity  plot  of  log  amplitude  and  the  absorption  coefficient  using  a  2D 
mesh  with  9664  nodes  for  a  source  and  detector  which  are  diagonally  opposite  to  each  other. 
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Table- 1  shows  the  RMS  error  with  respect  to  the  high  resolution  mesh  in  the  horizontal  cross- 
section  (as  indicated  by  the  dotted  line  in  Fig.  2)  using  the  method  described  earlier.  The  RMS 
error  calculated  here  was  also  calculated  along  different  cross-sections  of  the  model  and  a 
similar  trend  was  seen.  The  mesh  with  1785  nodes  was  found  to  have  an  RMS  error  of  less 
then  5%  as  compared  to  the  finest  mesh. 


Fig.  4.  Comparison  of  Jacobian  cross-section  with  respect  to  measurement  number,  (a).  The 
horizontal  cross-sectional  plot  of  the  sum  of  2D  circular  Jacobian  matrix  values,  from  center  to 
the  boundary  at  y  =  0  mm.  (b)  The  normalized  sum  of  2D  circular  Jacobian  matrix  values,  with 
respect  to  the  value  at  the  center  (at  x  =  0  mm,  y  =  0  mm).  The  legend  gives  number  of 
measurements  associated  with  each  plot. 

The  2D  mesh  with  1785  nodes  was  used  for  the  calculation  of  the  Jacobian  and  the 
expected  noise  level  in  the  amplitude  measurements  was  assumed  to  be  1%.  For  both  the 
heterogeneous  and  homogeneous  2D  cases,  the  number  of  useful  singular  values  above  the 
noise  level  were  calculated,  and  the  results  are  shown  in  Fig.  3(a).  Figure  3(b)  is  a  bar  chart 
showing  useful  measurements  in  percentage  [given  by  Eq.  (5)]  for  each  set  of  measurements. 
Figure  3(c)  is  a  plot  of  the  rank  of  the  Jacobian  versus  the  total  number  of  measurements  for 
meshes  having  different  resolution  starting  from  150  to  3569  nodes  versus  number  of 
measurements.  The  Jacobian  calculated  is  also  mapped  onto  a  reconstruction  (pixel)  basis 
ranging  from  10  x  10  to  40  x  40.  The  number  of  useful  singular  values  as  function  of  pixel 
basis  elements,  for  each  set  of  measurements,  are  plotted  in  Fig.  3(d).  Finally,  for  the  2D  case, 
Fig.  4  shows  the  total  sensitivity  distribution  at  the  mid-axis  cross-section,  as  a  function  of  the 
number  of  measurements.  Table  2  shows  the  number  of  useful  singular  values  of  the  3D 
model  Jacobian  which  are  above  the  noise  level  (1%)  for  the  three  different  strategies,  and 
indicates  the  effective  number  of  measurements  which  will  be  contributing  to  the 
reconstructed  image  space  and  quality.  The  number  of  useful  singular  values  is  higher  for  the 
three  layer  out-of -plane  strategy.  The  useful  percentage  of  measurements  is  higher  for  the  3D 
single  plane  of  data,  whereas  the  condition  number  is  very  high  for  the  3D  three-layer  out  of 
plane  case.  Similar  data  is  also  included  using  the  2D  circular  geometry  for  comparison 
purposes,  with  240  measurements  and  the  same  corresponding  optical  properties  as  the  3D 
model. 

The  plots  of  the  3D  Jacobian  magnitude  as  normalized  to  the  value  at  the  center  of  the 
model  are  shown  in  Fig.  5.  These  plots  shows  that,  all  the  three  strategies  of  data  collection  in 
3D  are  hypersensitive  (in  X  &  Y  direction)  at  the  boundary.  Moreover  this  is  pronounced  for 
the  3D  single -plane  case.  In  the  Z-direction  (not  shown)  it  was  found  that,  as  expected  that, 
the  sensitivity  decreases  as  the  position  moves  from  centre  to  boundary  for  all  the  three  cases. 
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Z=  0  mm  Z=  5  mm 


Z=  10  mm  Z=  15  mm 


Fig.  5.  The  normalized  cross-section  in  the  X-Y  plane,  showing  the  total  sensitivity  across  the 
dotted  line  in  Fig.  2,  from  x=  0  mm  to  x  =  43  mm  (center  to  boundary)  at  Y  =  0  mm  normalized 
with  respect  to  the  sensitivity  at  the  origin,  (i.e.  X  =  0,  Y  =  0&  Z  =  0  mm). 


Original  |ia  56  240  552  992 


1560  2256  3080  4032 


0.01  0.012  0.014  0.016  0.018  0.02 


Fig.  6.  The  reconstruction  of  the  pa  distribution,  using  noisy  simulated  data  of  log  amplitude, 
for  a  circular  object  with  an  absorbing  inhomogeneity  at  the  center.  Different  numbers  of 
measurements  were  used  as  denoted  above  each  image,  ranging  from  56  up  to  4032  data  points. 

The  forward  mesh  was  1785  nodes  and  the  pixel  basis  consisted  of  30x30  elements.  The 
original  pa  distribution  is  shown  as  the  first  image. 

The  2D  reconstruction  of  a  circular  object  with  a  centralized  absorption  anomaly  of 
diameter  of  10mm  using  different  number  of  measurements,  along  with  original  jj.a 
distribution,  is  shown  in  Fig.  6.  The  contrast  of  the  inhomogeneity  to  background  is  2:1  and 
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for  these  reconstructions  a  pixel  basis  of  30  x  30  elements  is  used,  with  a  forward  mesh 
consisting  of  1785  nodes.  Figure  7  shows  the  plot  of  logarithm  of  rms  error  in  the  horizontal 
cross-section  (as  sown  by  dotted  line  in  Fig.  6)  as  a  function  of  measurement  number.  The 
legend  of  the  figure  gives  the  position  of  the  inhomogeneity  (diameter  of  10mm). 

Table  3  summarizes  the  results  of  the  3D  reconstruction.  Figure  8(a)  shows  the 
reconstructed  absorption  coefficient  distributions  for  a  spherical  absorption  inhomogeneity 
(diameter  of  15  mm)  located  at  (0,  0,  0)  with  a  contrast  of  2:1  to  background,  using  the  data 
collected  from  the  three  strategies.  Figure  8(b)  shows  the  results  of  the  same  effort  with  a 
spherical  inhomogeneity  located  near  to  the  boundary  (30,  0,  0).  The  results  show  that  the 
quantitative  values  of  the  anomaly  increases  as  the  anomaly  is  moved  from  centre  to  boundary 
in  X  &  Y  direction.  The  anomaly  for  this  location  is  reconstructed  with  89%  quantitative 
accuracy  compared  to  the  15%  accuracy  for  central  location.  Finally  the  reconstructed 
absorption  coefficient  distribution  for  a  spherical  absorption  inhomogeneity  (diameter  -  15 
mm),  which  is  centered  at  (30,  0,  10)  are  shown  in  Fig.  8(c)  and  it  can  be  seen  that  single  layer 
case  reconstructed  the  anomaly  in  the  wrong  location.  Here,  both  the  in-plane  and  out-of- 
plane  strategies  are  able  to  give  up  to  84%  quantitative  accuracy  (Table  3). 


Fig.  7.  A  plot  of  logarithm  of  rms  error  in  the  horizontal  cross-section  of  |Ia  at  y  =  0  (as  shown 

in  original  |Ta  distribution  of  Fig.  6)  versus  number  of  measurements  for  various  positions  of  an 
absorption  inhomogeneity.  These  calculations  used  1785  nodes  in  the  mesh  of  the  forward 
problem  and  a  pixel  basis  of  30x30  elements  in  the  reconstruction. 

4  Discussion 

The  decrease  in  the  RMS  error  for  the  horizontal  cross-section  of  the  2D  Jacobian  for  a  given 
source -detector  (diagonally  opposite  each  other)  for  a  mesh  greater  than  1500  nodes  as 
compared  to  9664  nodes  (Table-1)  is  below  5%.  It  should  be  noted  that  the  other  kernels  of 
the  Jacobian,  for  example  J3  (30),  showed  better  accuracy  (2%)  when  the  mesh  had  1785 

dr 

nodes  or  greater.  As  with  many  iterative  reconstruction  problems,  optical  tomography  requires 
repeated  forward  calculations  and  re-computation  of  the  Jacobian,  thereby  increasing  mesh 
resolution  which  further  implies  increase  in  computational  time,  which  is  clearly  evident  from 
last  two  rows  of  Table  1.  A  computation  limit  of  10  seconds  per  iteration,  lead  to  a  choice  of 
mesh  resolution  with  1785  nodes  for  the  forward  problem  in  two-dimensional  case,  and 
extending  this  same  level  of  resolution  to  3D  would  require  nearly  80,000  nodes,  which  is 
near  the  limit  of  what  can  be  done  computationally.  Thus  much  of  the  2D  study  presented 
here  was  run  at  the  level  of  1785  nodes.  Since  the  computation  of  the  Jacobian  using  the  ITiM 
relies  on  the  discretization  of  the  domain  and  the  accuracy  of  the  numerical  model  depends  on 
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Table  3.  The  computation  time  and  accuracy  of  the  3D  reconstruction  is  shown  for  the  three 
different  data  collection  strategies,  along  with  three  different  locations  of  the  anomaly  for  each. 


Strategy 

Position  of 
anomaly 
(original) 

Iterati 

ons 

Total 

Computation 
time  (s) 

Quantitative 
accuracy  (%)  of 
the  reconstructed 
anomaly 

FWHM 
along  X- 
axis  (mm) 

FWHM 

along 

Z-axis 

(mm) 

3D:  1  layer 

(0,0,0) 

ii 

3179 

15% 

16.1 

25.2 

(30,0,0) 

14 

4046 

89% 

17.2 

23.3 

(30,0,10) 

10 

2890 

- 

- 

- 

3D  3layer  in- 
plane 

(0,0,0) 

14 

8022 

14% 

16.5 

25.3 

(30,0,0) 

14 

8022 

80% 

13.1 

18.7 

(30,0,10) 

12 

6876 

110% 

11.2 

18.6 

3D  3layer 
out-of-plane 

(0,0,0) 

6 

10926 

11% 

23.7 

24.1 

(30,0,0) 

9 

16389 

78% 

13.6 

18.9 

(30,0,10) 

8 

14568 

84% 

13.2 

18.7 

Fig.  8.  The  reconstructed  absorption  coefficient  distribution  for  the  cylindrical  object  with  a 
spherical  absorption  inhomogeneity  (diameter  of  15mm  and  contrast  2:1  with  respect  to 
background)  located  at  x,  y  and  z  locations  (a)  (0,0,0),  (b)  (30,0,0)  and  (c)  (30,0,10).  The  three 
columns  of  images  show  the  results  achieved  with  the  three  different  data  collection  schemes. 
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this  discretization  and  the  associated  integration  of  the  shape  functions,  the  resolution  of  the 
mesh  and  the  associated  optical  properties  will  affect  these  results.  For  example,  if  the 
absorption  coefficient  is  much  smaller,  then  lower  resolution  meshes  may  be  adequate,  as  the 
problem  becomes  more  energy  conserving,  whereas  for  a  higher  absorption  or  scattering 
problem,  a  higher  resolution  mesh  will  be  needed  to  ensure  numerical  accuracy  within  each 
FEM  element  for  a  lossy  problem.  Note  also  that  for  spectral  reconstruction  [3]  with  six 
wavelengths  data,  each  iteration  takes  about  30  sec  for  1785  nodes  mesh. 

For  a  heterogeneous  or  homogeneous  2D  case,  number  of  useful  singular  values,  which 
are  above  the  noise  level  (1%  in  amplitude)  showed  similar  trends  and  behavior  with 
increasing  numbers  of  measurements,  as  evident  from  Fig.  3(a).  Further,  the  percentage  of 
useful  measurements  (useful  singular  values)  drops  exponentially  as  the  number  of 
measurements  is  increased.  Fig.  3(b).  It  is  worth  noting  that  for  a  heterogeneous  model,  since 
light  propagation  becomes  more  complex,  and  in  this  case  more  diffusive,  the  total  number  of 
useful  measurements  is  slightly  lower  than  that  of  homogeneous  model.  In  this  work,  useful 
singular  values  are  defined  as  the  ones  which  are  above  noise  level  (1%).  This  is  used  only  for 
optimizing  the  parameters  used  in  the  reconstruction  procedure,  but  in  the  actual 
reconstruction  procedure,  regularization  is  used  to  reduce  the  condition  number. 

Next,  the  effect  of  the  2D  mesh  resolution  was  investigated,  for  it’s  impact  upon  the 
number  of  independent  available  measurements.  From  Fig.  3(c),  it  is  evident  that  if  the 
degrees  of  freedom  (mesh  resolution)  in  the  forward  problem  is  less  than  the  total  number  of 
measurements,  then  increasing  the  number  of  measurements  does  not  increase  the  number  of 
independent  measurements  (i.e.  the  rank),  since  the  rank  is  predominantly  restricted  by  the 
number  of  nodes  in  the  mesh.  For  example,  given  a  system  from  which  only  240 
measurements  are  available,  any  mesh  which  has  a  resolution  of  240  nodes  or  more  will  give 
the  same  number  of  independent  measurements.  Therefore  no  additional  measurements  can  be 
gained  in  terms  of  independent  information  by  increasing  the  mesh  resolution.  Given  a  2D 
mesh  of  1785  nodes,  for  example,  no  considerable  gain  in  independent  data  can  be  obtained 
when  the  number  of  measurements  are  increased  above  1560  (40  source  and  detectors).  At 
this  point,  it  will  be  worth  remembering  that,  in  real  time  there  is  a  physical  constraint  on 
number  of  measurements,  because  of  the  physical  geometry  and  fiber  size.  To  take  an 
example,  for  a  circular  test  phantom  of  86  mm  diameter  and  fiber  of  6  mm  diameter,  no  more 
than  40  fibers  (which  corresponds  to  1560  measurements)  can  be  arranged  around  the  outer 
boundary  of  domain.  However  this  issue  becomes  more  important  perhaps  for  non-contact 
imaging  systems  in  which  the  number  of  source-detector  locations  can  be  arbitrarily  large. 

Using  a  2D  mesh  of  1785  nodes,  the  effect  of  an  increase  in  the  reconstruction  (regular 
pixel)  basis  resolution  upon  reconstruction  was  investigated  [Fig.  3(d)].  An  increase  in  pixel 
basis  elements  increases  the  number  of  useful  singular  values,  but  there  is  no  significant 
improvement  in  the  pixel  basis  from  30x30  (900  elements)  to  40x40  (1600  elements).  This  is 
very  interesting,  since  one  would  assume  that  fewer  degrees  of  freedom  for  the  inverse 
problem  would  produce  a  better  solution.  But  although  the  problem  may  become  better  posed, 
the  rank  will  be  similar  to  that  shown  in  Fig.  3(d).  However,  these  results  indicate  the  best 
possible  resolution  obtainable  is  by  using  the  40  x  40  pixel  basis  and  again  these  results  will 
be  dependent  on  the  physical  problem  dimension  and  level  of  complexity.  Figure  4  shows  that 
increasing  the  number  of  measurements  for  a  2D  model  increases  the  sensitivity  of  the 
problem,  as  evident  from  magnitude  plot  of  the  Jacobian  (calculated  from  1785  nodes  mesh). 
Also  shown  in  Fig.  4  is  a  normalized  plot,  relative  to  the  central  value,  and  indicates  that  for 
fewer  number  of  measurements,  the  sensitivity  is  maximal  near  the  boundary  and  lower  at  the 
center,  as  expected.  By  increasing  the  number  of  measurements,  eventually  the 
hypersensitivity  near  to  the  boundary  reduced  and  the  sensitivity  became  uniform  regardless 
of  distance  from  boundary.  Finally,  it  is  observed  that  increasing  the  number  of  measurements 
above  552  (24  sources  and  detectors)  did  not  result  in  any  further  improvement  in  the 
sensitivity  distribution. 

For  the  3D  model.  Table  2  shows  that  three  layers  of  out-of-plane  measurements  yields  a 
higher  number  of  useful  singular  values,  but  the  useful  percentage  of  the  total  measurements 
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was  below  15%.  An  increase  in  number  of  measurements  means  more  data  acquisition  time 
and  more  computation  time.  Non-linear  iterative  image  reconstruction  procedures  in  NIR 
imaging  use  repeated  calculation  of  the  forward  data.  Therefore  increasing  the  number  of 
sources  and  measurements  substantially  increases  the  computation  time.  In  comparing  the 
three  layer  in-plane  and  three  layer  out-of-plane  data  collection  strategies,  having  more  than 
three  times  the  measurements  in  the  latter  case  improves  the  number  of  useful  singular  values 
only  by  22%.  The  improvement  in  the  number  of  useful  singular  values  is  not  significant  if  the 
data  acquisition  time  is  considered  as  well  as  the  computation  time.  The  magnitude  of  the 
singular  values  indicates  the  importance  of  that  eigenvector  in  the  image  space,  which  is 
directly  related  to  reconstructed  image  contrast  that  can  be  achieved.  To  compare  the 
magnitude  of  the  largest  singular  value,  even  though  it  is  at  its  highest  for  the  three  layer  out- 
of-plane  strategy,  it  should  be  noted  that  only  3  of  the  singular  values  are  above  164 
(magnitude  of  largest  singular  value  of  3D  31ayer-in-plane),  indicating  that  there  would  not  be 
dramatic  differences  in  the  reconstructed  image  contrast  in  both  these  cases.  If  the  magnitude 
of  largest  singular  value  in  2D  and  3D  are  compared,  in  2D  the  magnitude  is  higher,  whereas 
the  number  of  useful  singular  values  are  lower  than  3D,  indicating  that  the  modes  that 
contribute  to  image  space  are  fewer  and  the  quality  of  the  reconstructed  image  in  2D  will  be 
lower  than  3D.  Even  though  magnitude  of  the  singular  values  dictate  the  contrast,  the  singular 
vectors  associated  with  it  will  tend  to  affect  the  reconstructed  image  quality.  The  magnitude  of 
the  largest  singular  value  in  the  3D  31ayer  cases  are  the  same  because  of  the  smoothness  of  the 
singular  vectors  in  the  case  of  3D  3  layer:  out-of-plane,  the  reconstructed  image  quality  is 
better  than  the  rest  cases  (Fig.  8).  The  FWHM  analysis  also  confirms  this. 

It  should  be  noted  that  there  is  always  a  trade-off  between  image  quality  and  computation 
time.  Therefore  having  out-of-plane  data  increases  the  image  resolution,  but  taking  into 
consideration  the  overall  computation  time,  this  improvement  is  perhaps  not  so  significant. 
The  computation  time  per  iteration  is  high  in  the  case  of  out-of-plane  data  (computation  time 
per  iteration:  2D  problem  -  70  sec;  single -layer  -  289  sec;  three  layer:  in-plane  -  573  sec; 
three  layer:  out-of-plane  -  1821  sec). 

Figure  5  indicates  that  for  the  3D  model  with  a  single  measurement  plane  case,  the  total 
sensitivity  is  higher  near  the  boundary,  as  compared  to  the  three  plane  data  case  and  by 
increasing  the  number  of  measurements  the  sensitivity  near  the  boundary  is  decreased.  The 
results  show  that  although  the  sensitivity  is  still  higher  at  the  boundary  with  three  planes  of 
data  acquired,  there  is  no  significant  difference  in  the  sensitivity  pattern  observed  between 
three  layer  in-plane  or  out-of-plane  strategies. 

Since  only  one  component  of  the  full  Jacobian  matrix,  J2  in  Eq.  (2),  has  been  examined 
here,  images  have  also  been  reconstructed  for  |4a  using  log  amplitude  data  for  a  2D  forward 
mesh  of  1785  nodes  and  a  reconstruction  basis  30  by  30  pixel  basis.  Noisy  simulated  data 
were  generated  for  various  radial  positions  of  the  absorption  inhomogeneity  with  a  contrast  of 
2,  relative  to  the  background  and  having  a  diameter  of  10  mm.  The  log  of  RMS  error  was 
calculated  as  the  difference  in  the  original  and  the  reconstructed  horizontal  cross-sections  of 
each  image  (Fig.  6)  as  a  function  of  number  of  measurements  and  these  were  plotted  in  Fig.  7. 
The  results  show  that,  as  evident  from  Fig.  7,  although  there  is  a  decrease  in  the  RMS  error  as 
the  number  of  measurements  is  increased,  the  improvement  in  the  reconstructed  images  is  not 
significant  for  measurements  greater  than  552  (corresponding  to  24  fibers).  However,  for  a 
central  anomaly,  the  RMS  error  continued  to  decrease  with  increasing  number  of 
measurements,  whereas  for  an  anomaly  near  the  boundary  the  RMS  error  does  not  improve 
more  than  0.5%  with  respect  to  552  measurements. 

To  study  the  effect  of  data  collection  strategies  on  the  3D  reconstructed  image,  the  FWHM 
(Full  Width  at  Half  Maximum)  of  the  peaks  for  all  the  reconstructed  cases  have  been 
calculated  and  compared.  Table  3.  As  the  inhomogeneity  moves  from  the  centre  towards  the 
boundary,  the  FWHM  reduces  for  both  of  the  three  layer  cases  and  it  remains  approximately 
the  same  for  the  single  layer  case.  For  example,  when  the  inhomogeneity  is  placed  at  (30,0,0), 
Fig.  8(b),  the  FWHM  (in  the  X-cross  section)  values  for  single  layer  is  17.2mm  and  for  the 
three-layers  in-plane  and  out-of-plane  strategies  is  13.1mm  and  13.6mm  respectively.  It  is 
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evident  from  the  reconstruction  examples  that  the  quantitative  values  of  the  inhomogeneities 
increase  as  the  object  moves  from  the  centre  to  boundary,  which  is  in  close  match  with 
Jacobian  analysis  above.  Reconstruction  of  absorption  using  single  layer  data,  is  not  accurate, 
in  a  case  where  the  anomaly  is  not  presented  in  the  imaging  plane,  such  a  case  results  are 
presented  in  Fig.  8(c).  In  this  case,  single -layer  reconstructed  image  shows  the  inhomogeneity 
at  a  false  position  (reconstructed:  (30,0,0);  actual:  (30,0,10)).  Most  of  the  3D  NIR  studies 
indicate  that,  the  quantitative  accuracy  of  the  images  will  be  poor  due  the  partial  volume 
effect  in  three  dimensions[13, 16,17]  and  these  quantification  can  be  greatly  improved  by  the 
use  of  more  sophisticated  regularization  and  the  addition  of  penalty  terms  into  Eq.  (3). 

5  Conclusions 

In  this  investigation,  the  mesh  resolution  and  numerical  accuracy  in  the  2D  and  3D  forward 
problems  were  examined,  using  specific  data-collection  geometries.  Several  choices  such  as 
domain  size,  optical  properties  and  anomaly  position  and  size  were  kept  fixed,  relative  to 
typical  breast  cancer  imaging  situations.  It  was  shown  that  increasing  the  number  of 
measurements  increases  the  total  amount  of  information  available,  and  these  specifically 
enhance  the  recovery  of  the  central  region  of  the  model,  regardless  of  dimensionality.  Further, 
by  increasing  the  number  of  measurements,  the  rank  of  the  problem  (i.e.  amount  of 
independent  useful  information)  may  not  increase  if  the  degrees  of  freedom  (i.e.  number  of 
nodes  in  the  mesh)  are  low.  Reconstruction  basis  plays  an  important  role  in  the  inverse 
problem  and  it  has  been  found  that  a  pixel  basis  of  30  x  30  is  optimal  for  a  typical  breast 
imaging  problem. 

More  specifically  for  a  3D  imaging  problem,  this  work  has  shown  the  benefits  and 
drawbacks  of  multi-plane  data  collection  as  well  as  the  use  of  in-plane  versus  out-of-plane 
data  measurements  strategies.  It  has  been  shown  that  the  use  of  single -plane  of  data  in  a  3D 
model  is  perhaps  adequate,  in  terms  of  image  quality,  computation  time  and  data  collection 
time,  if  the  anomaly  being  imaged  is  within  the  plane  of  measurements.  Flowever,  if  prior 
information  such  as  plane  of  interest  is  not  known,  it  has  been  shown  that  multi-plane  data  is 
crucial.  The  use  of  in-plane  and  out-of-plane  data  has  been  addressed  and  is  shown  that 
although  the  use  of  out-of-plane  data  provides  more  independent  and  useful  information  for 
image  reconstruction,  the  magnitude  of  this  additional  information  does  not  provide  enough 
advantages  worth  the  data  acquisition  and  image  computation  time. 

Finally  it  is  worth  noting  that  the  3D  study  has  been  limited  to  16  source/detection  fibers 
per  plane.  The  addition  of  more  measurement  fibers  and/or  investigation  of  a  different  image 
reconstruction  basis,  such  as  those  performed  for  the  2D  problem  can  be  easily  extended  for 
the  presented  3D  problem.  The  technique  and  analysis  described  here  can  be  used  as  a  tool  to 
improve  resolution  and  contrast,  given  prior  information  about  the  domain  being  imaged.  This 
specific  study  was  undertaken  to  better  understand  the  parameters  and  capabilities  of  existing 
breast  imaging  system  at  Dartmouth  and  to  focus  on  software  improvements  which  may 
increase  its  recovery  of  lesion  information. 
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ABSTRACT 


Near-Infrared  (NIR)  tomographic  image  reconstruction  is  a  non-linear, 
ill-posed  and  ill-conditioned  problem,  and  so  in  this  study,  different  ways 
of  penalizing  the  objective  function  with  structural  information  were  inves¬ 
tigated.  A  simple  framework  to  incorporate  structural  priors  is  presented, 
using  simple  weight  matrices  that  have  either  Laplacian  or  Helmholtz 
structures.  Using  both  MRI-derived  breast  geometry  and  phantom  data,  the 
reconstructed  results  show  that  structural  priors,  which  would  be  available 
from  multi-modality  imaging,  give  superior  results  compared  to  having  no 
structural  prior  information.  Quantification  of  the  optical  properties  of  the 
tissue  is  more  accurate  and  the  structure  of  the  tissue  is  also  conserved  by 
incorporation  of  the  structural  information.  More  over,  parameter  reduction 
(i.e.  hard  prior  information)  in  the  imaging  field  through  the  enforcement  of 
spatially  explicit  regions  gives  erroneous  results  (over-estimates  the  absorption 
coefficient  by  360%  and  under-estimates  the  scattering  coefficient  by  88%), 
if  the  structural  prior  information  from  one  of  the  regions  is  imperfect 
by  as  little  as  7%  of  the  area.  Even  with  the  most  accurate  priors,  this 
parameter-reduction  technique  can  over-estimate  the  scattering  coefficient  of 
the  tumor  by  over  100%  as  estimated  by  experimental  studies.  Incorporation 
of  less-restrictive  prior  (soft-prior)  information  can  be  achieved  with  either 
Helmholtz  or  Laplacian  structured  regularization  matrices.  Using  multi¬ 
layered  phantom  data,  it  is  shown  that  the  Helmholtz  type  of  regularization 
matrix  gives  the  best  estimate  of  the  scattering  parameter  and  the  Laplacian 
provides  best  estimate  for  the  absorption  parameter.  It  is  also  proven  that  for 
smaller  values  of  the  wave  number,  the  Helmholtz  and  Laplacian  structures 
give  similar  estimates  for  optical  properties.  In  the  same  framework,  it  is  also 
demonstrated  that  applying  unreasonable  constraints  to  the  imaging  problem 
amplify  the  noise  in  the  solution,  leading  into  physiologically  invalid  estimates 
of  optical  properties. 


2 


1.  Introduction 


Near  Infrared  (NIR)  optical  tomography  is  a  method  that  uses  light  in  the  650-900nm 
wavelength  range  to  recover  images  of  the  internal  spatial  distribution  of  tissue  optical 
properties,  absorption  (or  chromophore  concentrations)  and  scattering  parameters.1-3  When 
imaged  at  multiple  wavelengths,  these  quantities  can  be  used  to  estimate  tissue  hemoglobin 
and  water  concentration3  and  have  the  advantage  of  being  acquired  non-invasively  and 
without  ionizing  radiation.  The  imaging  procedure  can  be  rapidly  or  repeatedly  applied  to 
investigate  physiological  state,  and  systems  can  be  integrated  into  conventional  imaging 
platforms  such  as  X-ray  mammography,  Ultrasound  and  MRI.4-9  These  hybrid  systems 
have  been  shown  to  achieve  superior  performance  in  terms  of  resolution  and  quantitative 
accuracy  which  should  provide  more  accurate  physiological  data  from  the  tissue  under 
investigation.4-16  However,  a  fundamental  question  is  how  to  utilize  the  spatial  information 
from  the  clinical  system  optimally  to  maximize  the  accuracy  of  NIR  tomography.  In  this 
study,  the  ability  to  improve  the  quantitative  accuracy  of  regions  imaged  with  NIR  to¬ 
mography  was  investigated,  in  the  setting  where  prior  spatial  information  is  readily  available. 

This  work  explores  image  reconstruction  strategies  that  take  advantage  of  multimodality 
image  data,  in  particular  combination  of  MRI  with  NIR  optical  tomography  for  breast 
cancer  imaging.  MRI  provides  structural  information  at  high  spatial  resolution  (near 
1  mm),  whereas  NIR  imaging  has  relatively  poor  spatial  resolution  (near  4-7  mm). 
Yet  MR  imaging  would  benefit  from  the  molecular-specific  signatures  available  through 
NIR,3,16-18  specifically  tissue  hemoglobin  content,  oxygenation  level,  and  water  as  well  as 
scattering  particle  size  and  number  density.19  The  inverse  problem  (image  reconstruction 
procedure)  in  NIR  imaging  is  known  to  be  a  non-linear,  ill-posed  and  ill-conditioned 
problem.20  Use  of  structural  information  in  NIR  reconstruction  schemes  has  been  explored 
by  several  research  groups.  For  example,  Li  et.  al7  have  used  the  data  derived  from  X-ray 
mammography  for  choosing  different  regularization  parameters  for  the  region  of  interest 
(ROI)  and  surrounding  tissue,  and  have  shown  that  the  contrast  and  resolution  of  the 
reconstructed  images  can  be  improved.  Srinivasan  et  al.21  have  developed  a  three-step 
reconstruction  process  for  improving  the  quantification  accuracy  of  small-objects  in  NIR 
tomography,  where  they  use  the  conventional  NIR  reconstructed  images  (first  step)  as  a 
structural  prior  for  the  last  two  steps.  Earlier  papers  have  shown  that  optical  contrast 
can  be  correlated  to  MR  contrast6,9,13  and  that  structural  MRI  images  can  be  used  to 
reduce  the  number  of  unknown  parameters  to  be  estimated.14  The  difficulty  with  parameter 
reduction  approaches  (referred  to  as  hard-priors)  is  the  potential  of  introducing  error  by 
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imposing  incorrect  model  assumptions  and  introducing  variation  due  to  uncertainty  in 
the  prior  information  (even  when  the  underlying  model  is  appropriate).  For  example,  the 
features  which  lead  to  contrast  in  one  imaging  system  may  not  be  spatially  coregistered 
with  those  that  produce  contrast  in  another  imaging  system.  Further,  segmentation  of 
congruent  features  always  includes  classification  errors  resulting  from  digitization.  Recently, 
Boverman  et.  al10  have  shown  that  even  imperfect  priors  which  encode  breast  background 
structure  improves  anomaly  localization,  but  at  the  expense  of  biasing  the  spectroscopic 
dimension  of  the  image  reconstruction.  Another  type  of  approach,  often  described  as 
soft-priors,  for  constraining  the  problem,  penalizes  the  variation  within  regions  which  are 
assumed  to  have  the  same  properties  by  controlling  regularization.  Brooksby  et.  al15, 16,22 
have  developed  a  Laplacian  type  of  regularization  that  allows  intra-region  variability.  This 
is  a  method  which  works  well  even  if  the  confidence  in  the  prior  structural  information  is  low. 

This  paper  develops  a  more  generalized  framework  for  incorporating  the  structural  priors 
into  the  NIR  image  reconstruction  process  and  explores  a  covariance-based  constraint  scheme 
adopted  from  finite  differencing  of  the  Helmholtz  equation  in  addition  to  the  soft  and  hard 
prior  approaches  noted  above.  The  former  allows  optical  property  variation  with  a  given 
region,  reducing  biases  caused  by  the  use  of  imperfect  prior  information.  The  results  indicate 
that  imperfect  structural  information  can  generate  errors  in  the  hard-priors  case,  whereas 
the  soft-priors  are  able  to  quantify  regions  more  appropriately.  Simulation  and  experimental 
studies  are  performed  to  demonstrate  the  superior  reconstruction  image  quality.  These  types 
of  procedures  are  needed  to  improve  NIR  imaging  both  in  terms  of  high  spatial  resolution 
available  from  MRI  and  high  contrast  inherent  in  the  NIR  signal. 

2.  Methods 

2. A.  Diffusion-based  Light  Transport  Model 

Light  transport  in  breast  tissue  can  be  described  accurately  by  the  Diffusion  Equation  (DE), 
which  is  an  approximation  to  the  Radiative  Transport  Equation  (RTE).23  In  the  frequency- 
domain,  the  DE  is  given  by 

— V.k(r)V<f>(r,u;)  +  (/xa(  r)  +uj/c)$(t,uj)  =  Q0(r,u;)  (1) 

where  <3>(r,u;)  is  the  photon  density  at  position  r  and  light  modulation  frequency  is  given  by 
c o  (in  this  work,  c a  =  100  MHz).  The  isotropic  source  term  is  represented  by  Q0(r,u)  and 
speed  of  light  in  tissue  by  c.  /xa(r)  is  the  optical  absorption  coefficient  and  k(r)  is  the  optical 
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diffusion  coefficient,  which  is  defined  as 

k(r)  =  — - — — — - — -  (2) 

3k(r)  +  /i'(r)] 

where  /4(r)  is  the  reduced  scattering  coefficient,  which  is  defined  as  fj,'s  =  /xs(l  —  g).  ps  is 
the  scattering  coefficient  and  g  is  the  anisotropy  factor.  A  Robin-type  (Type-III)  boundary 
condition  is  applied  to  exactly  model  the  refractive- index  mismatch  at  the  boundary.24  The 
boundary  data  for  a  frequency  domain  system  are  the  amplitude  and  phase  of  the  measured 
signal,  which  is  used  with  a  Finite  Element  Method  (FEM)  based  reconstruction  procedure 
to  obtain  the  internal  spatial  distributions  of  ga  and  /i( . 

2.B.  Standard  image  reconstruction 

The  objective  function  (If)  for  this  procedure  can  be  written  as 

min 

0=  ,  {||y-F(pa,k)||2  +  A||(pa,k)-(pao,ko)||2}  (3) 

Ra,  k 

Where,  ’F’  is  the  forward  operator  that  generates  the  model  response  and  ’y’  is  the  exper¬ 
imental  measured  data.  ||.||2  represents  L2-Norm  of  the  vector.  This  is  also  known  as  the 
Tikhonov  approach,26  where  ’A’  is  the  regularization  parameter  that  balances  the  degree  to 
which  differences  between  the  current  estimate  of  optical  properties  and  the  starting  values 
are  allowed  and  data- model  misfit.  More  specifically,  it  is  the  ratio  of  the  variances  of  the 
data  noise  and  the  parameter  (A  =  o2y/crf^a  kj)-  Minimizing  Eq.  3  (by  setting  first  derivatives 
with  respect  to  pa  and  k  to  zero)  leads  to  an  update  equation 

( JTJ  +  A ,J)(5//a,  5k)  =  JT(y  -  F(pa,  k))  -  A[(/za,  k)  -  (pa0,  k0)]  (4) 

where,  ’J’  is  the  Jacobian  matrix  and  I  is  the  identity  matrix.  Note  that  JT  J  is  ill-conditioned; 
however  I  stabilizes  the  matrix.  However,  a  slight  deviation  from  this  update  equation  is 
generally  employed,  which  is  also  known  as  Levenberg-Marquardt  (LM)  type  of  regularization 
procedure,27,28  assuming  [(S/ia,8k)  =  (fia,  k)  —  (pa0,k0)]  leading  to29 

(■ JTJ  +  2A/)(5/ia,  5k)  =  JT( y  -  F(pa,  k))  (5) 

Most  of  the  literature  reports  A*  =  2A, 20,29  which  is  true  only  for  Levenberg-Marquardt 
minimization,  which  does  not  involve  the  parameter  field  in  the  objective  function  (Eq. 
3)  ,2'’28  In  this  LM  approach,  A*  typically  starts  being  the  ratio  of  the  variances  and  is 
reduced  at  each  of  the  iterations  by  a  small  factor  (in  here,  it  is  y/lO  and  also  multiplied 
by  the  maximum  of  the  diagonal  values  of  JTJ).  The  iterative  procedure  is  repeated  until 
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experimental  data  matches  with  modeled  data  within  a  preset  value  e  (~  data  noise  level). 
In  general,  the  initial  values,  (yUao,k0),  are  obtained  from  a  pre-reconstruction  step  where 
the  data  is  calibrated  by  a  homogeneous  fitting  procedure.30,31 

2.  C.  Inclusion  of  a  priori  information 

The  objective  function  with  inclusion  of  prior  information  is  given  as15 

min 

0=  ,  { lly  F(/xa,  k)  || 2  +  A||L[(/ra,  k)  —  (/ra0,  k0)]  ||2}  (6) 

ha,  k 

Here  also  A  is  the  ratio  of  the  variance  of  the  data  noise  to  parameter  held  and  L  is  a  penalty 
matrix  (dimensionless  in  all  the  cases  considered  in  this  work)  which  can  be  derived  from 
MRI  structural  priors  as  indicated  below.  The  update  equation  resulting  from  this  procedure 
is: 

(JTJ  +  XLTL)(6iaa,Sk)  =  JT{ y  -  F(/za,k))  -  ALTL[(/ra,  k)  -  (//. a0,k0)]  (7) 

In  this  work,  each  location  in  the  computational  discretized  model  is  labeled  according  to 
tissue  type  (fatty,  hbroglandular  or  tumor)  determined  from  MRI  Tl-weighted  images.15,16,22 
It  is  also  assumed  that  there  is  no  covariance  between  the  different  regions  of  the  imaging 
domain.  Since  the  domain  model  does  not  itself  change  throughout  the  iterative  reconstruc¬ 
tion  algorithm,  the  L-matrix  is  calculated  before  the  reconstruction  procedure  and  it  is  used 
through  out  the  process  to  penalize  the  solution.  Two  forms  for  the  L-matrix  considered 
in  this  work  are  discussed  in  the  subsections  below,  including  the  Laplacian  and  Helmholtz 
structures. 


2.C.I.  a  Laplacian  Structured  Regularization  Matrix 
The  Laplace  equation  in  parameter  u(r)  can  be  written  as 

V2u(r)  =  0  (8) 


A  finite  difference  approximation  to  this  equation  for  ’N’  number  of  equi-space  (step  size  = 
h)  nodes  can  be  written  as32 

V2w(r)/?2  «  ui  +  u2  +  ■  ■ .  —  Nun/2  +  . . .  +  uN-i  +  uN  =  0  (9) 


Dividing  the  whole  equation  ’-N’  leads  to 


V2w(r) 


-u  i 


N 


+ 


-U2 

N 


+  .  .  .  +  1ijv/2  +  • 


—UN-1  —Un 

'  +  N  +  ~W~ 


=  0 


(10) 
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The  L  matrix  is  a  matrix  that  relates  each  nodal  property  of  the  numerical  model  to  all  other 
nodes.  Therefore  given  a  node  i  within  the  mesh,  its  relationship  to  another  node  j  having 
Laplacian  structure  (Eq.  10)  within  the  same  mesh  can  be  given  as15,16 

!0  if  i  and  j  are  not  in  the  same  region 

—  1/N  if  i  and  j  are  in  the  same  region  (11) 

1  if  i  =  j 

where  N  is  the  number  of  finite  element  mesh  nodes  comprising  a  given  region.  In  this 
case,  LTL  approximates  a  second-order  Laplacian  smoothing  operator  within  each  region, 
and  functionally  works  to  average  the  update  within  a  region,  while  allowing  discontinuity 
between  different  regions. 


2.C.2.  a  Helmholtz  Structured  Regularization  Matrix 

The  Helmholtz  equation  in  parameter  u(r)  for  a  damped  wave  can  be  written  as 

V2w(r)  —  n2u{r )  =  0  (12) 


where  k  is  the  wave  number,  specifically  k,  =  1//,  where  l  covariance  length.32  I  also  represents 
the  decay  length  scale  over  which  the  parameter  u(r)  has  correlation.  Making  k  =  0,  will 
give  the  Laplace  equation  (Eq.  8).  A  finite  difference  approximation  to  this  equation  for  ’N’ 
number  of  equi-space  (step  size  =  h)  nodes  can  be  written  as32 

(V2  -  k2)  u(r)h 2  «  Mi  +  u2  +  ■  ■  ■  +  [—(IV  +  ( Kh)2)]uN/2  +  . . .  +  uN- 1  +  uN  =  0  (13) 

Dividing  the  whole  equation  by  ’  —  (N  +  (k/i)2)’  gives 


(V2  -  k2)  u(r)  : 


—  Ml _ _ -M2 

-(V  +  (k/i)2)  -(N  +  (Kh.)2) 


-  UN/  2  ■ 


— MJV-1 


— UN 


-( N  +  (nh )2)  -(N+(kIi)2) 


=  0 


(14) 


Writing  this  as  a  generalized  L-matrix  form  similar  to  Eq.  11 


0  if  i  and  j  are  not  in  the  same  region 

L(i,j )  =  <  —  JV+(1h.it)2  if  i  and  j  are  in  the  same  region  (15) 

|  1  if  i  =  j 

For  the  FEM  nodes  case,  h  is  chosen  to  be  the  distance  between  the  nodes.  More  over, 
k  —  l/l  is  generally  chosen  to  be  the  inverse  of  the  size  of  the  feature  (tumor  in  this  case) 
in  the  imaging  domain.  This  case  is  also  known  as  Best  Prior  Estimate  (BPE).  As  the  prior 
structural  information  is  available  through  MRI,  /  is  chosen  to  be  the  diameter  of  the  target 
(tumor)  in  the  BPE  case.  In  this  case,  LTL  (L  given  by  Eq.  15)  approximates  a  second-order 
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Helmholtz  smoothing  operator.  To  determine  the  effect  of  k  on  the  parameter  reconstruction, 
different  values  for  k  are  chosen.  It  is  shown  that  for  small  values  of  k,  which  corresponds 
to  a  large  correlation  length  (/),  both  Laplacian  and  Helmholtz  structures  recover  the  same 
optical  property  distribution. 

2.D.  Breast  geometry  -  Effect  of  imperfect  a  priori  information 

The  techniques  described  in  Sec.  2  were  used  to  reconstruct  images  from  synthetic  data 
with  1%  random  noise  added.  Numerical  experiments  using  synthetic  data  generated  on  a 
volunteer  MRI  Tl-weighted  breast  images  with  incorrect  priors  to  show  the  effectiveness  of 
soft-priors.  Figure  1(a)  (first  column)  shows  the  original  distribution  of  three  tissue  layers, 
namely  fatty  (/ia  =  0.006  mm-1  and  fi's  =  0.6  mm-1),  fibro  glandular  (/za  =  0.012  mm-1  and 
p!s  =  1.2  mm-1)  and  tumor  (/ia  =  0.018  mm-1  and  fi's  =  1.8  mm-1)  for  the  breast  geometry 
(labeled  as  0,  1  and  2  respectively,  in  Fig.  1(a)  first  column).  Sixteen  light  collcction/delivery 
fibers  were  arranged  equally  spaced  on  a  circle  (indentions  in  Fig.  1(a)).  In  succession,  one 
fiber  was  used  as  the  source  while  all  other  fibers  served  as  detectors  which  provided  a 
total  of  240  measurements.  In  these  studies,  the  source  was  modeled  as  a  Gaussian  profile 
with  a  Full  Width  Half  Maximum  (FWHM)  of  3  mm  to  most  accurately  represent  the 
light  applied  using  the  clinical  system  used,  and  was  placed  at  a  depth  of  one  transport 
scattering  distance  from  the  tissue  boundary.  A  mesh  of  1785  nodes  (corresponding  to  3418 
linear  triangular  elements)  was  used  for  the  diffusion  model  and  reconstruction  calcula¬ 
tions.33  A  total  of  7%  of  the  glandular  layer  (label-1)  FEM  nodes  were  labeled  (relative 
to  the  original  glandular  layer  nodes)  as  fat  (label-0)  to  introduce  imperfect  structural  priors. 

The  same  initial  estimates  (optical  properties  of  region  ’O’)  were  used  as  homogeneous 
starting  conditions.  The  iterative  procedure  was  stopped  once  the  data-modcl  misfit  (resid¬ 
ual)  did  not  improve  by  more  than  2%  when  compared  with  the  previous  iteration.  The 
starting  value  for  A  is  chosen  to  be  25000  and  75  for  pa  and  k  respectively,  derived  from  the 
noise  characteristics,  for  Eq.  5.  The  same  values  are  chosen  for  Eq.  7. 

2.E.  Breast  geometry  -  Effect  of  data  noise  level 

The  techniques  described  in  Sec.  2  were  used  to  reconstruct  images  from  synthetic  data  with 
1,  3,  5  and  10%  Gaussian  noise  to  see  the  effect  of  data  noise  level  on  the  reconstruction 
techniques.  The  breast  geometry  (and  the  optical  properties)  were  equivalent  to  the  previ¬ 
ous  section  (Fig.  1(a)),  however,  perfect  spatial  priors  were  used.  The  same  FEM  mesh  as 
described  above  was  employed  in  the  forward  and  reconstruction  problems.  Optical  proper¬ 
ties  of  region  ’O’  were  used  as  initial  guess  for  the  iterative  procedure.  The  regularization 


parameter  (A)  and  stopping  criterion  was  chosen  according  to  each  data  noise  level. 

2.F.  Phantom  studies 

A  multi-layered  gelatin  phantom  (86  mm  diameter,  25  mm  height)  was  fabricated  with 
different  optical  properties  using  heated  mixtures  of  water  (80%),  gelatin  (20%)  (G2625, 
Sigma  Inc.),  India  Ink  for  absorption  and  TiCU  (titanium  oxide  powder,  Sigma  Inc.)  for 
scattering15  were  solidified  by  cooling  to  room  temperature  (see  Fig.  3).  Different  layers 
of  gelatin  were  constructed  by  successively  hardening  gel  solutions  containing  different 
amounts  of  ink  and  TiCU-  A  cylindrical  hole  (diameter  of  16  mm  and  height  of  24  mm) 
was  filled  with  liquid  to  mimic  the  tumor.  The  first  column  in  Fig.  4  shows  the  axial 
cross-section  of  three-layers  of  the  phantom  (Fig.  3)  where  the  region  labeled  ’O’  has  the 
homogeneous  optical  properties  (/za  =  0.0065  mm-1  and  /i's  =  0.  65  mm-1),  similar  to 
the  typical  fatty  layer  in  the  breast6  and  a  thickness  of  10  mm.  The  hbroglandular  layer 
(diameter  76nnn)  also  has  homogeneous  optical  properties  (region  labeled  T’)  of  fia  =  0.01 
mm-1  and  /%  =  1.0  mm-1.  The  tumor  (represented  by  the  region  labeled  ’2’)  has  a  diameter 
of  16  mm  with  optical  properties  of  [ia  =  0.02  mm-1  and  /i(  =  1.2  mm-1.  The  optical 
properties  were  validated  by  measuring  large  cylindrical  samples  of  each  layer.  Appropriate 
mixtures  of  Intralipid  and  India  ink  were  used  to  achieve  the  desired  optical  parameters  of 
the  tumor.  Data  was  acquired  using  a  clinical  NIR  system34  where  the  fibers  were  marked 
and  photographed  to  extract  region  information  (analogous  to  MRI  images).  This  regional 
information  was  used  to  label  the  corresponding  regions  in  the  FEM  mesh.15  A  mesh  of  1785 
nodes  (corresponding  to  3418  linear  triangular  elements)  was  used  for  the  diffusion  model 
calculations  and  a  mesh  having  1360  nodes  was  used  in  the  reconstruction.33  NIR  data 
was  calibrated  using  a  reference  homogenous  phantom  to  obtain  initial  optical  properties 
estimates  and  minimize  the  variation  between  the  16  optical  channels  according  to  standard 
practice  in  human  imaging  studies.30,31 


3.  Results 

Reconstructed  /ia  and  /%  images  obtained  from  the  noisy  simulated  data  with  imperfect 
(7%  error)  glandular  layer  priors  using  the  methods  described  in  Sec.  2  are  shown  in  Fig.  1. 
Using  hard  priors,  the  total  number  of  unknowns  is  reduced  to  6  parameters  (/aa  and  /%  for 
each  of  the  3  regions)  and  these  images  are  presented  in  the  last  column  of  Fig.  1(a).  The 
images  from  two  different  approaches  of  soft-priors  are  shown  in  the  middle  2  columns;  the 
Erst  column  shows  the  expected  results.  For  the  Helmholtz  case,  av  =  1/8  mm  (BPE)  was 
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used,  where  8  mm  is  the  diameter  of  the  tumor.  Cross-sectional  plots  of  the  reconstructed 
fia  and  /it  distributions  along  the  dotted  line  in  Fig.  1(a)  (first  column)  are  provided  in 
Fig.  1(b).  Table  1(a)  and  (b)  show  the  mean  and  standard  deviation  of  the  optical  property 
estimations  in  each  region  of  the  reconstructed  images.  Note  that  the  NIR  reconstruction 
procedure  without  prior  information  (not  shown)  did  not  generate  meaningful  images  in 
this  complex  case. 

Figure  2(a)  shows  the  reconstructed  images  from  the  data  with  5%  noise  using  all  four 
techniques  described  in  Sec.  2.  The  first  column  of  Fig.  1(a)  shows  the  actual  distribution 
of  optical  properties.  Figure  2(b)  shows  the  mean  and  standard  deviation  values  (as  error 
bars)  of  reconstructed  images  using  different  techniques  for  different  regions  of  the  breast 
geometry  with  increasing  data  noise  level.  In  the  Helmholtz  case,  k  —  1/8  mm  (BPE) 
was  used.  The  actual  values  are  also  plotted  for  the  comparison.  It  can  be  clearly  seen 
Laplacian  regularization  gives  lesser  standard  deviation  for  the  absorption  coefficient 
(/ia)  reconstruction  compared  to  the  Helmholtz  structure.  On  the  other  hand,  Helmholtz 
structure  performs  better  than  Laplacian  in  the  case  of  scattering  coefficient  (/i(). 

A  photograph  of  the  phantom  used  to  collect  data  at  785nm  with  16  fibers  in  a  single  plane 
(giving  240  measurements)  is  shown  in  Fig.  3.  Images  obtained  from  the  procedures  described 
in  Sec.  2  are  presented  in  Fig.  4(a)  along  with  cross-sectional  plots  of  the  reconstructed  /ja 
and  ii's  distributions  in  Fig.  4(b).  Table  2(a)  and  (b)  show  the  mean  and  standard  deviation 
of  the  optical  property  estimations  in  each  region  of  the  reconstructed  images.  Here  also,  for 
Helmholtz  case,  the  BPE  (k  =  1/16  mm)  was  used.  Figure  5(a)  gives  the  results  for  different 
values  of  n  (given  on  the  top  of  the  each  column,  true  distribution  is  shown  as  first  column  in 
Fig.  4(a))  in  the  Helmholtz  case.  Corresponding  cross-sections  are  plotted  in  Fig.  5(b).  The 
mean  and  standard  deviation  values  for  each  of  this  case  are  also  given  in  Table  2(a)  and  (b). 


4.  Discussion 

The  reconstructed  results  (Fig.  1,  2,  4  and  5)  show  that  the  structural  priors  improve  the 
reconstructed  image  quality  dramatically.  The  penalized  problem  formulation  (different 
type  of  regularizations)  generates  smoother  images  resulting  in  smaller  standard  deviations 
from  the  mean  values  (see  Table-2)  as  compared  to  the  generalized  problem  that  does  not 
incorporate  prior  information. 
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The  hard-prior  case  produces  significant  optical  property  value  error  when  the  structural 
a  priori  information  is  imperfect  in  the  breast  geometry  (Fig.  1  and  Table-1).  In  this  case, 
a  7%  variation  in  the  definition  of  the  glandular  layer  caused  false  estimates  of  the  optical 
properties.  On  the  other  hand,  soft-priors  (Laplacian  and  Helmholtz)  yield  good  estimates 
for  each  layer.  Hard-priors  over-estimate  the  tumor  absorption  coefficient  by  360%  and 
under-estimate  its  scattering  coefficient  by  88%  (Fig.  1).  Soft-priors  are  within  the  6%  of 
the  expected  values  even  with  the  error  in  the  structural  prior.  Note  that  in  this  particular 
case,  hard-priors  failed  when  7%  of  glandular  layer  made  as  a  fatty  layer,  yet  below  this 
error  value  it  gave  reasonable  estimates  of  optical  properties. 

With  perfect  structural  priors,  increasing  the  data  noise  level  also  increases  the  quan¬ 
tification  error  in  the  reconstructed  images  (Fig.  2(a)  and(b)).  Hard-Priors  clearly  fail 
in  quantifying  the  tumor  optical  properties  as  the  data  noise  level  increases.  Soft-Priors 
does  better  in  the  quantification  than  Hard-Priors.  It  is  also  evident  that  incorporation 
of  structural  information  is  key  for  accurate  quantification  of  the  optical  properties.  The 
experimental  results  (Fig.  4)  from  the  three-layer  gelatin  phantom  (Fig.  3)  also  show  that 
incorporation  of  perfect  structural  information  improves  the  quantification  and  quality  of 
the  reconstructed  images.  Specifically,  the  mean  and  standard  deviation  of  the  reconstructed 
optical  property  values  in  each  region  are  both  more  accurate  and  more  precise  where  the 
priors  are  included.  The  mean  values  (Table-2)  show  that  the  absorption  coefficient  for  the 
Region-0  (fat)  is  under-estimated  by  77%  in  the  case  of  the  Helmholtz  regularization  (BPE). 
This  error  can  be  explained  by  the  fact  that  the  Euclidean  distance  between  the  nodes 
was  used  rather  than  the  distance  between  the  nodes  along  the  boundary  of  a  particular 
region.  Both  Laplacian  and  Hard-Priors  over  estimate  the  scattering  coefficient  of  the 
tumor  region  by  a  factor  of  2.  It  is  known  that  the  photon  path  length  is  affected  by  the 
scattering  coefficient.  By  constraining  the  problem  based  on  the  distance,  one  can  expect  to 
estimate  the  scattering  better.  More  over,  the  Helmholtz  equation  allows  wave  propagation, 
which  models  the  photon  diffusion  better.  Table- 2  indicates  that  the  Helmholtz  technique 
always  produces  more  quantitative  accuracy  for  the  scattering  coefficient  estimation  and  the 
Laplacian  technique  is  best  for  the  absorption  coefficient  estimation  (as  well  as  Fig.  2  and  4). 

Theoretically  Helmholtz  and  Laplacian  structures  are  identical  when  k,  =  0  (equivalently  l 
is  large).  Figures  5(a)  and  (b)  (as  well  as  Tablc-2)  show  that  when  n  =  1/86  mm  (/  =  86  mm 
is  the  diameter  of  the  domain),  Laplacian  and  Helmholtz  structures  give  reasonably  close 
reconstruction  values  of  optical  properties,  which  indicates  the  expected  trend  presented  in 
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this  paper  for  the  two  methods.  It  also  indicates  that  the  BPE  case  (k  =  1/16  mm)  gives 
the  best  results,  as  the  priors  applied  are  close  to  the  true  structure  of  the  feature  (tumor). 
These  results  also  indicate  that  unreasonable  constraints  (like  k  —  1/5  mm,  Fig.  5(a)  first 
column)  makes  the  estimation  problem  amplify  the  noise  resulting  in  physiologically  invalid 
(scattering  coefficient  is  greater  for  fatty  layer  compared  to  fibroglandular  layer)  estimates 
of  optical  properties.  In  the  tumor  region,  as  hi  decreases,  the  estimated  values  of  optical 
parameters  become  closer  to  expected  values  (Table- 2  and  Fig.  5,  hi  =  1/43  and  hi  =  1/86 
cases).  This  is  due  to  the  correlation  length  becoming  larger,  making  the  covariance  in  the 
neighboring  nodes  larger. 

In  this  study,  it  is  shown  that  imperfect  priors  (commonly  caused  by  improper  image 
segmentation  and  image  artifacts  in  MRI  or  X-ray)  can  lead  to  error-prone  results  in  the 
hard-prior  case  whereas  soft-priors  are  more  immune  to  uncertainty  in  the  prior  data.  It  is 
also  shown  that  the  techniques  used  to  incorporate  the  soft  structural  prior  information  in¬ 
fluences  the  image  outcome,  which  may  lead  to  improvements  in  image  accuracy  if  properly 
implemented.  Srinivasan  et.  al35  have  found  that  5%  error  in  optical  properties  introduces 
approximately  45%  error  in  the  chromophore  concentration  when  a  limited  number  of  wave¬ 
lengths  are  used  for  imaging.  The  correct  ’’soft”  inclusion  of  a  priori  information  therefore 
can  be  expected  to  lead  to  a  more  accurate  quantification  of  chromophore  concentrations  as 
well.  It  should  be  noted  that  over-weighting  of  the  penalty  term  in  the  problem  formulation 
may  make  the  solution  ignore  the  data-model  misfit  and  emphasize  smooth  feature  extrac¬ 
tion.  The  techniques  developed  in  this  work  were  applied  for  two-dimensional  test  objects, 
and  can  be  easily  extended  to  three-dimensional  case.  A  more  extensive  study  of  this  is  left 
for  future  investigations. 

5.  Conclusions 

This  work  has  investigated  several  ways  of  incorporating  structural  information  into  an  itera¬ 
tive  image  reconstruction.  The  results  have  been  supported  by  gelatin  phantom  experiments 
that  represent  multi-layered  structure  which  is  commonly  found  in  breast  tissue  with  adipose 
(fatty)  tissue  on  the  exterior  and  fibroglandular  tissue  nearer  to  the  interior.  Soft-priors 
allow  the  tissue  optical  properties  to  vary  within  predefined  regions,  unlike  hard-priors  which 
constrain  these  zones  to  be  homogeneous.  Hard-priors  were  found  to  perform  poorly  when 
the  prior  information  contained  area  errors  as  small  as  7%  which  can  easily  be  produced 
by  most  segmentation  algorithms.  True  boundary  extraction  from  MRI  images  intro¬ 
duces  unavoidable  segmentation  and  discretization  errors  that  are  better  tolerated  when  the 
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structural  information  is  encoded  through  the  soft-prior  approach  involving  a  penalty  matrix. 


The  results  reported  here  indicate  that  the  optical  properties  of  different  tissue  types  can 
be  quantified  more  accurately  when  their  estimation  is  properly  guided  by  ’’soft”  structural  a 
priori  information.  The  problem  formulation  and  results  presented  in  this  work  indicate  that 
data  from  other  imaging  modalities  such  as  ultrasound  or  x-ray  tomography,  could  also  be 
used  as  the  source  of  the  structural  prior.  In  the  cases  investigated,  the  Helmholtz  structure 
always  gives  a  better  estimation  of  scattering  coefficient.  However,  the  Laplacian  type  of 
regularization  leads  to  more  superior  absorption  coefficient  estimate.  So  it  is  reasonable  to 
conclude  that  Laplacian  structure  gives  the  best  estimates  of  total  hemoglobin  concentration 
(. Hbr ),  hemoglobin  oxygen  saturation  (StO-2%)  and  water  fraction  ( H-20 )  (which  are  the 
main  absorbers).  Helmholtz  structure  gives  the  best  estimates  of  the  scattering  power  and 
scattering  amplitude  (scattering  parameters).  The  framework  presented  here  can  also  be 
extended  to  other  regularization  terms  such  as  total  variation  minimization  or  spectral  prior 
constraints,  which  may  be  studied  in  future  work. 
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List  of  Tables 


Table  1:  Mean  and  standard  deviation  of  the  reconstructed  (a).  /ia  and  (b).  h's  values  in 
different  regions  (labeled  in  first  column  of  Fig.  1(a))  recovered  with  simulated  data  having 
1%  random  noise  and  imperfect  structural  information  defining  region  ’1’  (7%  reduction 
compared  to  the  original  segmentation).  The  corresponding  reconstructed  images  are  shown 
in  Fig.  1(a) 

Table  2:  Mean  and  standard  deviation  of  the  reconstructed  (a).  fia  and  (b).  h's  values 
in  different  regions  (labeled  in  first  column  of  Fig.  4(a))  recovered  from  the  experimental 
phantom  data.  The  corresponding  reconstructed  images  are  shown  in  Fig.  4(a)  and  5(a). 

List  of  Figures 

Figure  1:  (a).  Actual  /ja  and  fi's  distributions  from  a  breast  (obtained  from  a  volun¬ 
teer)  are  shown  in  the  first  column.  Optical  properties  for  the  region  labeled  ’O’  (fat)  are: 
Ha  =  0.006  mm-1  and  Hs  —  0.6  mm-1.  Region  ’1’  (hbroglandular)  values  are:  Ha  =  0.012 
mm-1  and  //(,  =  1.2  mm-1.  Region  ’2’  (tumor)  values  are:  Ha  —  0.018  mm-1  and  h's  = 
1.8  mm-1.  Reconstructed  Ha  and  h's  images  from  different  techniques  with  simulated  data 
having  1%  random  noise  and  imperfect  structural  information  defining  region  ’1’  (7% 
reduction  compared  to  the  original  segmentation)  are  shown  in  the  rest  of  the  columns. 
The  middle  two  columns  use  soft  prior  structural  information  while  the  last  row  shows  the 
result  with  hard  prior  information.  In  the  Helmholtz  case,  k  —  1/8  mm  (BPE)  was  used, 
(b).  Cross-sectional  plots  along  the  dotted  line  in  the  actual  image  (see  first  column  of  (a)) 
of  the  true  and  reconstructed  Ha  and  h's  distributions. 

Figure  2:  (a).  Reconstructed  Ha  and  h's  images  from  different  techniques  with  simu¬ 
lated  data  having  5%  random  noise  and  perfect  structural  priors  (actual  images  are  shown 
in  the  Erst  column  of  Fig.  1(a)).  The  first  column  shows  the  reconstruction  results  without 
the  use  of  prior  information.  The  middle  two  columns  use  soft  prior  structural  information 
while  the  last  row  shows  the  result  with  hard  prior  information.  In  the  Helmholtz  case,  k  = 
1/8  mm  (BPE)  was  used.  (b).  The  mean  values  and  standard  deviations  (plotted  as  error 
bars)  in  Ha  and  h's  f°r  different  regions  of  breast  geometry  (labeled  in  actual  image)  with 
increasing  noise  level  (1%  to  10  %). 
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Figure  3:  Photograph  for  gelatin  phantom  (representing  the  idealized  two-dimensional 
cross-sectional  geometry  shown  as  first  column  in  Fig.  4(a))  used  in  the  experimental  studies. 

Figure  4:  (a).  Actual  /ia  and  /i's  distributions  (axial  cross-section)  of  phantom  (Fig. 
3)  case  are  shown  in  the  first  column.  Optical  properties  for  the  region  labeled  ’0’  are:  fia  = 
0.0065  mm-1  and  /i(  =  0.  65  mm-1.  Region  ’1’  values  are:  / ia  =  0.01  mm-1  and  fi's  =  1.0 
mm-1.  Region  ’2’  (tumor)  values  are:  /ia  =  0.02  mm-1  and  /i(  =  1.2  mm-1.  Reconstructed 
Ha  and  n's  distribution  from  different  techniques  (discussed  in  Sec.  2)  from  the  experimental 
phantom  data.  Second  column  of  images  does  not  use  prior  information.  The  middle  rows 
use  soft  prior  structural  information  and  the  last  row  of  images  were  recovered  with  hard 
priors.  In  the  Helmholtz  case,  n  =  1/16  mm  (BPE)  was  used.  (b).  Cross-sectional  plots  along 
the  dotted  line  in  the  actual  image  (see  first  column  of  (a))  of  the  true  and  reconstructed 
/ia  and  fi!s  distributions. 

Figure  5:  (a).  Reconstructed  /ia  and  //s  images  from  the  experimental  phantom  data 
using  Helmholtz-type  regularization  matrix  for  different  values  of  k,  which  are  given  at  the 
top  of  each  column,  (b).  Cross-sectional  plots  along  the  dotted  line  of  the  actual  images  in 
Fig.  4(a)  (first  column)  are  shown  with  the  data  from  reconstructed  fia  and  /i(  images  in 
(a).  The  best  prior  estimate  (BPE)  case  (k  =  1/16  mm)  is  also  presented  for  comparison. 
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Tables 


Methods 

Region- 0 

Region- 1 

Region-2 

Actual 

0.006 

0.012 

0.018 

Laplacian 

0.0064±0.0010 

0.0117±0.0018 

0.0156±0.0018 

Helmholtz  ( k  =  1/8) 

0.0062±0.0011 

0.012±0.002 

0.0156±0.0017 

Hard  Priors 

0.006 

0.0118 

0.0843 

(a) 


Methods 

Region- 0 

Region- 1 

Region-2 

Actual 

0.6 

1.0 

1.8 

Laplacian 

0.63±0.10 

1.13±0.18 

1.67±0.23 

Helmholtz  (k  —  1/8) 

0.64±0.09 

1.09±0.16 

1.64±0.22 

Hard  Priors 

0.64 

1.13 

0.23 

(b) 


Table  1 
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Methods 

Region- 0 

Region- 1 

Region-2 

Actual 

0.0065 

0.01 

0.02 

No  Priors 

0.0022±0.0005 

0.0061±0.0032 

0.0192±0.0044 

Laplacian 

0.0031±0.0002 

0.0051±0.0005 

0.0174±0.0029 

Helmholtz  (k  =  1/16) 

0.0015±0.0005 

0.0058±0.0009 

0.0241±0.0043 

Hard  Priors 

0.0032 

0.005 

0.0213 

Helmholtz  (k  =  1/5) 

0.0009±0.0006 

0.0061±0.0008 

0.0191±0.0031 

Helmholtz  (k,  =  1/43) 

0.0027±0.0003 

0.0052±0.0007 

0.0234±0.0043 

Helmholtz  (k  =  1/86) 

0.0022±0.0005 

0.0061±0.0032 

0.0192±0.0044 

(a) 


Methods 

Region- 0 

Region- 1 

Region-2 

Actual 

0.65 

1.0 

1.2 

No  Priors 

0.64A0.40 

0.66A0.27 

0.76A0.16 

Laplacian 

0.38A0.03 

0.63A0.07 

2.37A0.41 

Helmholtz  (k  =  1/16) 

0.46A0.02 

0.59±0.02 

1.08A0.12 

Hard  Priors 

0.37 

0.63 

2.74 

Helmholtz  (k  =  1/5) 

0.60A0.01 

0.57A0.01 

0.82A0.06 

Helmholtz  (k  —  1/43) 

0.39±0.03 

0.63A0.03 

1.19±0.13 

Helmholtz  (k  —  1/86) 

0.39±0.03 

0.63A0.04 

1.21±0.14 

(b) 


Table  2 
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Figures 


Figure  1(a) 
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Figure  1(b) 
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Figure  2(a) 
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ABSTRACT 


Diffuse  Optical  Tomography  (DOT)  involves  estimation  of  tissue  optical  prop¬ 
erties  using  noninvasive  boundary  measurements.  The  image  reconstruction 
procedure  is  a  non-linear,  ill-posed,  and  ill-determined  problem,  so  overcoming 
these  difficulties  requires  regularization  of  the  solution.  While  the  methods 
developed  for  solving  the  DOT  image  reconstruction  procedure  have  a  long 
history,  there  is  less  direct  evidence  on  the  optimal  regularization  methods, 
or  exploring  a  common  theoretical  framework  for  techniques  which  uses 
least-squares  (LS)  minimization.  A  generalized  least-squares  (GLS)  method 
is  discussed  here,  which  takes  into  account  the  variances  and  covariances 
among  the  individual  data  points  and  optical  properties  in  the  image  into 
a  structured  weight-matrix.  It  is  shown  that  most  of  the  least-squares 
techniques  applied  in  DOT  can  be  considered  as  special  cases  of  this  more 
generalized  LS  approach.  The  performance  of  three  minimization  techniques 
using  the  same  implementation  scheme  is  compared  using  test-problems  with 
increasing  noise  level,  and  increasing  complexity  within  the  imaging  field. 
Techniques  that  use  spatial-prior  information  as  constraints  can  be  also 
incorporated  into  the  GLS  formalism.  It  is  also  illustrated  that  inclusion 
of  spatial  priors  reduces  the  image  error  by  at  least  a  factor  of  2.  The 
improvement  of  GLS  minimization  is  even  more  apparent  when  the  noise  level 
in  the  data  is  high  (as  high  as  10%),  indicating  that  the  benefits  of  this  ap¬ 
proach  are  important  for  reconstruction  of  data  in  a  routine  setting  where  the 
data  variance  can  be  known  based  upon  the  signal  to  noise  of  the  instruments. 
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1.  Introduction 


Image  reconstruction  methods  used  in  diffuse  optical  tomography  (DOT)  are  mainly 
dependent  on  the  type  of  data,  diffuse  light  model  and  amount  anatomical/spectral 
priors  available.  There  are  numerous  reconstruction  techniques  available  in  the  literature 
depending  on  the  application.1-5  Yet  despite  the  volume  of  work  in  this  area  there  is  no 
single  investigation  with  a  direct  comparison  of  the  LS  minimization  techniques  using  the 
same  implementation  scheme,  especially  in  terms  data  noise  level  and  complexity  in  the 
test  holds.  Most  of  the  comparisons  in  the  literature  have  been  in  terms  of  implementation 
of  minimization  and  convergence  rates  of  one  or  two  techniques  at  hand.1-5  This  work 
addresses  this  problem  and  compares  minimization  methods  (more  specifically  different 
types  of  regularization)  with  the  same  implementation  scheme  for  a  direct  quantitative 
comparison.  Moreover,  usage  of  weight  matrices  in  the  regularization  which  include  the 
variance  and  covariance  properties  of  data  and  image  space  is  extensively  explored  here.  A 
new  covariance  form  borrowed  from  meteorological  studies  is  introduced  and  proven  to  be 
effective  for  reconstructing  highly  noisy  data  in  the  generalized  theoretical  frame  work. 

Near  Infrared  (NIR)  DOT  involves  reconstructing  images  of  optical  properties  from 
transmission  measurements  using  wavelengths  from  650nm  to  lOOOnm  to  interrogate 
tissue.1,6-8  Optical  absorption  and  scattering  images  obtained  using  multiple  wavelengths 
can  be  used  to  estimate  tissue  hemoglobin,  water  concentration  and  scattering  amplitude 
and  scattering  power.8  To  overcome  the  inherent  low-spatial  resolution  in  DOT,  there  is 
a  considerable  interest  in  developing  hybrid  systems,9-27  which  use  the  spatial  mapping 
of  one  system  as  the  template  for  DOT.  Image  formation  from  the  data  collected  by 
these  (stand-alone/hybrid)  systems  involves  solving  an  inversion  problem.  This  paper 
describes  least-squares  (LS)  minimization  techniques  to  solve  the  inverse  problem,  and  to 
quantitatively  compare  their  performance  in  a  systematic  series  of  simulations.  The  inverse 
problem  (image  reconstruction  procedure)  in  DOT  is  known  to  be  a  non-linear,  ill-posed, 
and  ill-determined  problem,2  and  to  solve  such  a  problem,  a  regularization  term  must  be 
added  to  constrain  the  solution  space  in  order  to  obtain  a  meaningful  image.  There  are 
many  regularization  methods  available  in  the  literature,  and  this  work  focuses  on  the  fact 
that  most  LS  techniques  presented  in  the  literature  can  be  encompassed  within  a  generalized 
theoretical  framework,  which  includes  a  regularization  matrix  that  is  based  upon  weights 
from  the  data  and  parameter  variances.  Note  that  appendix-A.l  gives  the  terminology  used 
in  this  work  along  with  definitions  of  symbols. 
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Because  of  the  interest  in  using  spatial  information  derived  from  conventional  imaging 
modalities  in  the  DOT  inverse  problem,  a  number  of  methods  have  been  presented  in  the 
literature.9-27  These  techniques  were  initially  proposed  by  Barbour  et  al9and  Schweiger  et 
al13  and  used  in  to  improve  the  quantitative  outcome  of  reconstructed  images.  Ntziachristos 
et  al14  used  the  MR  information  to  divide  the  imaging  domain  into  tumor  and  non-tumor 
regions  to  make  the  problem  better  posed.  Li  et  al1'  used  an  X-ray  tomosynthesis  volume  to 
segment  the  breast  into  different  sub-regions  and  used  different  regularization  parameters 
depending  on  the  size  of  the  sub-regions.  Recently  Guven  et  al24  proposed  a  Bayesian  frame 
work  to  include  spatial  prior  information  in  an  effective  way  which  will  not  bias  the  image 
reconstruction  problem  to  imperfect  anatomical  priors.  Pogue  and  Paulsen,10  Brooksby 
et  al,18, 21,25  Yalavarthy  et  al26  have  extended  these  approaches  for  the  use  of  anatomical 
prior  information  in  which,  depending  on  the  connectivity  and  size  of  the  sub-region,  the 
regularization  term  was  scaled.  Even  though  the  effect  of  imperfect  spatial  prior  information 
on  the  image  reconstruction  is  a  very  active  research  area,23, 24,26  it  was  assumed  here  that 
the  spatial  priors  were  perfect.  Other  ongoing  studies  are  examining  this  more  complex  issue. 


2.  DOT  Forward  Problem 

DOT  involves  solving  a  model  (forward)  and  estimation  (inverse)  problem,  sequentially 
as  illustrated  in  Fig.  1.  In  this  section,  the  forward  problem  is  described,  which  involves 
generating  the  measurement  data,  for  a  given  set  of  optical  property  estimates  within  the 
tissue,  using  a  finite  element  solution  to  the  diffuse  transport  equation. 


Light  propagation  in  a  turbid  elastic-scattering  media,  like  tissue,  is  treated  as  ‘neutral- 
particle  transport’  rather  than  ‘wave-propagation’  and  in  the  frequency-domain,  the  Diffusion 
Equation  (DE)  is  used,  which  is  given  by2,28 

-  V.-D(r)  V<L(r,  a;)  +  (na( r)  +  iu/c)$(r,  c o)  =  Q0( r,  uj)  (1) 


where  <f>(r,u;)  is  the  photon  density  at  position  r  and  the  light  modulation  frequency  is  given 
by  u  (ui  =  27rf,  in  this  work  f  =  100  MHz).  The  isotropic  source  term  is  represented  by 
Q0(r,cu)  and  the  speed  of  light  in  tissue  by  c,  which  is  constant  here.  /xa(r)  is  the  optical 
absorption  coefficient  and  D{ r)  is  the  optical  diffusion  coefficient,  which  is  defined  as 


D(  r)  = 


(2) 


3[/ia(r)  +/4(r)] 

where  /4(r)  is  the  reduced  scattering  coefficient,  which  is  defined  as  /i's  =  /xs(l  —  g).  /is  is  the 
scattering  coefficient  and  g  is  the  anisotropy  factor.  A  Robin  (Type-III)  boundary  condition 
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is  applied  to  model  the  refractive-index  mismatch  at  the  boundary.29  The  measured  data 
for  a  frequency  domain  system  are  the  amplitude  and  phase  of  the  transmitted  signal.  If 
F  is  the  forward  model  (Finite  Element  Method  (FEM)  in  here)  which  gives  the  fluence  at 
every  point,  then  the  modeled  data  G{y)  can  be  obtained  by  sampling  the  forward  model  at 
the  boundary  given  internal  spatial  distributions  of  optical  properties  and  source-detector 
locations,  where  /j  represents  the  parameters  to  be  estimated  (ji—  [-D(r);  /ia(r)]). 

G(/x)  =  S{F(y)}  (3) 

The  details  of  the  FEM  formulation  of  the  forward  model  are  given  in  Refs. 30-32  The 
results  presented  are  restricted  to  frequency- domain  data,  more  specifically  data  (y)  is 
the  natural  logarithm  of  the  amplitude  (A)  and  phase  ( 9 )  of  the  frequency-domain  sig¬ 
nal.  Defining  A  and  9  in  terms  of  modeled  data,  A  =  ^ J Re{G(y )}2  +  Lm{G(p)}2  and 
9  =  tan ~4  (Im{G(n)}/ Re{G(/a)}).  The  Jacobian  (J),  which  gives  the  rate  of  change  of 
modeled  data  with  respect  to  parameters,  is  calculated  using  the  adjoint  method.30  Even 
though  the  actual  parameters  being  estimated  are  D(r)  and  /aa(r),  the  results  are  presented 
in  terms  of  /ia(r)  and  y's( r),  which  are  spectroscopically  more  meaningful. 

3.  Least-Squares  Minimization  Techniques 

This  section  outlines  several  different  minimization  schemes  used  in  this  work.  These 
techniques  are  used  to  solve  the  inverse  problem  (Fig.  1),  which  is  achieved  by  minimizing 
the  objective  function  (D)  over  the  range  of  p.  Minimizing  the  objective  function  can  be 
achieved  by  several  different  approaches.  The  most  common  approaches  involve  obtaining 
repeated  solutions  of  the  forward  model  and  recomputation  of  the  Jacobian  (J)  (and  its 
inversion)  at  every  iteration  because  of  the  non-linear  nature  of  the  problem.  There  are 
also  gradient-based  optimization  schemes  available  in  the  literature33,34  to  minimize  the 
objective  function  which  does  not  require  an  explicit  inversion  of  the  Hessian  matrix.  In 
here  direct  methods,  known  as  full-Newton  approaches,2  are  employed  in  minimization  for 
all  the  regularization  techniques  used  for  a  fair  comparison.  LS  minimization  has  the  effect 
of  reducing  high  frequency  noise,  leading  to  smooth  images  of  optical  properties.  Total 
variation  methods  and  variants  of  this  are  used  to  obtain  edge  preservation  in  reconstructed 
images.27,35  Solving  the  inverse  problem  using  LS  minimization  can  be  also  seen  from 
Bayesian  prospective  to  get  maximum  a  posteriori  (MAP)  estimate.24,36,37  A  correlation 
between  the  Bayesian  frame  work  and  LS  minimization  techniques  is  given  in  Refs.,12,38,39 
but  usage  of  Bayesian  frame  work  requires  one  to  choose  a  particular  noise  model  for 
both  data  and  image  space,  which  might  not  reflect  the  actual  noise  characteristics  unless 
some  prior  information  is  available.  Here,  the  focus  is  on  Least-Squares  (LS)  minimization 
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techniques  with  a  focus  on  what  the  value  of  the  regularization  method  can  be.  The  LS 
methods  are  divided  into  two  groups;  (1).  Without  spatial  priors  (2).  With  spatial  priors. 

Without  Spatial  Priors 

3. A.  Levenberg-Marquardt  (LM)  Minimization 

This  approach  is  also  known  as  a  trust-region  method5,39  where  experimental  data  is  matched 
with  modeled  data  iteratively.40,41  The  objective  function  for  the  DOT  problem  is  defined  as 

tt  =  {||y-G'Gu)||2}  (4) 

where  y  is  the  data  and  G(p)  is  the  modeled  data.  This  equation  is  minimized  by  setting 
the  first-order  derivative  equal  to  zero. 


First-Order  Condition :  Minimizing  0  with  respect  to  p,  which  is  achieved  by  setting 

f  =  0 

0[l 


on 

dp 


Jt5  =  0 


(5) 


where  6  is  the  data-model  misfit,  6  =  y  —  G(p),  J  is  the  Jacobian,  T  represents  the  matrix 
transpose  operator. 


Iterative  Update  Equation :  Imagine  a  sequence  of  approximations  to  p  represented  by 
Pi,  then  using  Taylor  series  on  G(pi )  and  expanding  around  /p_i  gives 

G(pi)  =  G(pi-i)  +  J Api  +  . . .  (6) 


where  A p^  =  p^  —  pi-\.  Rewriting  5  utilizing  the  first  two  terms  of  Eq.  6  (ignoring  the  rest, 
equivalently  linearizing  the  problem)  gives 


5i  =  y  -  G(pi )  =  y  -  G(pi_ i)  -  JAp{  =  5 -  JApi 


Rewriting  Eq.  5  for  the  ith  iteration 
Substituting  Eq.  7  into  Eq.  8  gives 


JT5i  =  0 


JT(di- 1  -  JAp^  =  0 


(7) 

(8) 

(9) 


Further  simplification  leads  to  the  update  equation: 


A  pi 


(10) 
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When  ,JTJ  is  ill-conditioned,  a  diagonal  term  is  added  to  stabilize  the  problem.  In  this  case, 
the  update  equation  becomes: 


J  -) -  oi  I 


A/q 


(11) 


Where  A/q  is  the  update  for  the  parameter  in  the  ith  step.  Note  that  a  monotonically 
decreases  with  iterations  (always  >  0),  and  also  that  a  >  ||5||2.  The  iterative  method  (or  its 
modified  version)  is  the  commonly  used  minimization  technique  in  DOT.  It  can  be  seen  from 
Eq.  10  and  11,  when  a  becomes  zero  in  Eq.  11  it  becomes  Eq.  10.  It  is  also  important  to  note 
that  JT  J  is  always  symmetric,  because  (JT J)1  =  -JT(JT)T  =  -JT  J.  The  advantage  of  using 
this  method  is  in  the  simple  choice  of  a  regularization  parameter  (cc).  The  limitations 41  of 
this  method  include: 


•  ,JT  J  must  be  positive  definite. 


•  The  initial  guess  (/u,0)  should  be  close  to  the  actual  solution. 

•  The  update  equation  (Eq.  11)  does  not  solve  the  first-order  conditions  unless  a  =  0. 

•  Since  parameters  are  not  involved  in  the  minimization  scheme,  the  inverse  problem 
may  be  unstable. 

Even  though  J 1 J  is  not  positive  definite  in  DOT,  the  LM  approach  (or  its  modified  version) 
has  been  used  successfully  in  a  number  of  instances.2, 6,7,28,42 


3.B.  Tikhonov  Minimization 

The  generalized  objective  function43,44  in  the  Tikhonov  case  includes  parameters  in  the 
minimization  function,  which  is  defined  as: 

^  =  {||y-G(/i)||2  +  A||L(q-q0)ir}  (12) 

where  A  is  the  Tikhonov  regularization  parameter  and  L  is  a  dimensionless  regularization 
matrix  (in  this  work).  Here,  /i 0  is  the  prior  estimate  of  the  optical  properties,  which  in  DOT 
has  typically  been  obtained  from  calibrating  the  data.45,46 


Choice  of  A:  Rewriting  Eq.  12,  normalizing  both  terms  by  their  variances  yields 

0  =  Jy  -  g(r)II2  ,  \\l{t-^)\\2 

1  K)2  o)2 

where  ay  is  the  standard  deviation  in  the  data  y  and  cr/4_/i0  is  the  standard  deviation  in 
the  optical  properties  (or  deviation  from  the  prior  estimate  of  optical  properties).  Note  that 
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the  variance  of  data-model  misfit  (6  =  y  —  G(n ))  is  assumed  from  the  data,  i.e.  (o^)2  = 
(<xy)2  +  with  (oG(/i))2  =  0  because  synthetic  data  was  used.  Multiplying  Eq.  13  by 

Gy2  and  comparing  the  result  with  Eq.  12  leads  to: 


A  = 


by, 


o-,/- 


(14) 


which  shows  that  the  Tikhonov  regularization  parameter  (A)  should  be  equal  to  the  square 
of  the  ratio  of  the  standard  deviation  in  data  to  the  standard  deviation  of  the  parameters. 
This  is  a  subtle  yet  important  point,  especially  since  this  parameter  is  rarely  defined  this 
way,  and  is  most  commonly  derived  empirically. 


First-Order  Condition-.  Minimizing  fl  with  respect  to  /i,  which  is  achieved  by  setting 

on 


f2  =0 

Ofl 


^  =  jTs  -  \LTL(fi  -  no)  =  o. 

<7/1 


(15) 


Update  Equation :  Rewriting  Eq.  15  for  the  \th  iteration  leads  to 

JT$i  -  \LTL(ni  -  Ho)  =  0  (16) 

Substituting  Eq.  7  into  Eq.  16  results  in 

JT(di- 1  —  J  A  Hi)  —  \LT  L(Hi-i  +  Ahz  —  Ho)  =  0.  (17) 


Further  simplification  leads  to  the  iterative  update  equation: 


J 1 J  +  AL 1 L 


Z\/i;  —  J  <5j— i  AL  L(/li_i  Ho)' 


(18) 


Note  that  LTL  is  symmetric.  The  constraint  on  the  choice  of  L  is  that  it  must  be  positive 
definite.44  In  the  absence  of  spatial  priors,  a  common  choice  for  the  form  of  L  is  the  identity 
matrix  (I),  which  leads  to  the  update  equation 


J  + AI 


A  Hi 


JTh;_ 


-Mhi-l  —  ho) 


(19) 


Refer  to  appendix-A.2  for  an  analysis  of  the  Tikhonov  regularization  in  terms  of  singular 
values.  This  regularization  method  is  particularly  common  for  ill-posed  problems.  The  ad¬ 
vantage  of  the  method,  is  that  it  includes  parameters  within  the  minimization  scheme  which 
can  be  selected  to  improve  the  stability  of  the  solution.  Its  limitations  are  that 


•  it  requires  a  prior  opinion  about  the  noise  characteristics  of  the  parameter  and  data 
spaces  (for  A),  and 


•  it  does  not  take  into  account  the  individual  variances  of  the  data  points/parameters, 
nor  their  covariances. 

However,  the  simplicity  of  the  approach  makes  it  attractive  for  use  in  ill-posed  problems. 
When  the  dynamic  range  of  the  data  is  large  (as  in  DOT),  incorporation  of  the  maximum 
variance  in  the  data  will  cause  the  minimization  to  bias  the  solution  to  specific  data  points 
(e.g.  near  the  boundaries  at  source-detector  locations  in  DOT).  To  reduce  the  effect  of  bias, 
one  can  employ  a  generalized  least  squares  (GLS)  minimization  scheme,  described  in  the 
next  section. 

3.C.  Generalized  Least  Squares  (GLS)  Minimization 

Generalized  least  squares  minimization  schemes  have  been  proposed  in  the  context  of 
Tikhonov  minimization  in  the  literature,1,5,38  in  which  there  is  some  ambiguity  in  choos¬ 
ing  the  regularization  parameter  (A).  In  here,  a  direct  inclusion  of  weight  matrices  (which 
are  inverses  of  covariance  matrices)  in  the  minimization  scheme  was  employed  to  explic¬ 
itly  remove  the  dependence  of  reconstructed  image  quality  on  the  choice  of  regularization 
parameter.  This  type  of  choice  leads  to  an  objective  function:47,48 

n  =  {(y  -  G(p,))TWs(y  -  G(/i))  +  (/i  -  iiofW^iv  -  ii0)}  (20) 

where  Ws  is  the  weight  matrix  for  data-model  misfit  (5)  with  Ws  =  ( cov(6 ))~4  (Appendix 
A-4  of  Ref.47).  f'T/t-/J0  is  the  weight  matrix  for  optical  properties  (/ u-juo )  with  I'TM-/i0  = 
( cov(n  —  /i0))_1  (Appendix  A-4  of  Ref.47).  Explicit  forms  for  these  weight  matrices  are 
discussed  below.  Since  both  are  inverses  of  covariance  matrices,  they  are  symmetric  and 
positive  definite. 


First-Order  Condition :  Minimizing  D  with  respect  to  /i,  which  is  achieved  by  setting 
=  0  produces 

—  =  JTW56  -  -  ho)  =  0.  (21) 


Update  Equation :  Similar  to  Tikhonov  approach,  linearizing  the  problem  leads  to  the  iterative 
update  equation:48 


JtW5J  +  W^_M0  =  JTW5hi_i  -  -  fi0) 


(22) 
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3.C.I.  Choice  of  W$ 


Since  simulated  data  was  used  here,  in  the  formation  of  the  weight  matrix  (covariance  ma¬ 
trix),  it  was  assumed  that  the  cov(S)  is  due  to  measurement  error  only,  which  yields47 

Ws  =  [cou(<S)]_1  =  [cov(y  -  G(^))]-1  =  [< cov(y )]~7  (23) 

where  cov  represents  the  covariance  operator.  In  the  simulation,  typically  one  generates  the 
forward  data  and  adds  noise  to  it  to  form  synthetic  data. 


y  =  G()i)  +  Oyt J 


(24) 


Where  r/  is  a  random  number  vector.  Typically,  a  random  number  generator  which  follows 
a  normal  distribution  with  zero  mean  and  unity  variance  is  used.  Here,  ay  is  the  standard 
deviation  of  the  data,  assuming  the  noise  is  totally  uncorrelated  (white  noise)  in  which  case, 
the  covariance  matrix  becomes47 


cov{y)itj 


0  if  i  ±  j 

( °y)i  if  i  =  j 


(25) 


Since  synthetic  data  was  used  in  this  paper,  the  weight  matrix  for  the  data  {W$)  becomes 
diagonal.  In  the  experimental  case,  one  needs  to  collect  an  ensemble  of  data  sets  from  which 
a  covariance  matrix  can  be  computed.  In  this  case,  ’N’  data  sets  needs  to  be  collected  using 
the  same  phantom  (different  homogeneous  phantoms  need  to  be  used  for  different  signal 
levels),  where  ’N’  needs  to  be  a  large  number.  From  this  ensemble  of  {y} 


{y}  =  y  +  {y}  (26) 

where  y  is  the  true  or  mean  value  of  data  and  {y}  is  perturbation  due  to  noise.  This  leads 
to 

_  iiiiT 

[«*>(</)]  =  [®»(i/)]  =  ^  (27) 

substituting  27  in  Eq.  23  gives  Wg.  Note  that  in  Eq.  23,  it  was  assumed  the  cov(S)  is  clue 
to  measurement  error,  which  is  also  true  in  the  case  of  experimental  data,  as  the  data  is 
calibrated  to  remove  the  offset  and  match  the  modeled  data.45,46 


3.C.2.  Choices  of  W^,-^ 

Here,  two  forms  were  considered  to  highlight  the  versatility  of  the  procedure,  even  though 
many  other  forms  of  the  covariance  matrix  can  exist. 
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Analytical  Covariance  (AC)  form :  Borrowed  from  the  meteorological  studies,  assuming 
the  parameter  held  obeys  the  Helmholtz  equation,  an  analytical  form  (for  one- dimensional 
infinite  domain  case)  for  the  covariance  matrix  is:47 

[cov(n  -  p0)]  =  K-w,)2  (X  +  y)  e”^  (28) 

where  rt]  is  the  separation  distance  between  locations  and  /  is  the  correlation  length  scale. 
(cqt_M0)2  is  the  expected  variance  for  p  —  p0-  In  this  case,  the  weight  matrix  is  constructed 
from  Wy^  =  {cov{n  -  p0))_1 


Local  Laplacian  (LL)  form :  Here,  the  weight  matrix  is  formed  directly  using  a  local 
Laplacian  operator5, 49,50  between  neighboring  locations,  where  Hy/J0  =  (1  /(aM_w)2)MTM, 
where  M  is  the  local  Laplacian  matrix,  which  is  defined  as 


Mij 


0 

-1 

IE. j  Mij 


if  i  and  j  are  not  neighbors 
if  i  and  j  are  neighbors 

if  i  =  j 


(29) 


Computation  of  Hy,^  requires  an  estimate  of  variance  of  parameters  ((a^.y2),  same  is 
true  for  calculation  of  Tikhonov  regularization  parameter  (Eq.  14).  The  expected  variance 
can  be  computed  in  many  ways,  the  most  common  way  for  imaging  problems  is  from  the 
literature  (which  gives  the  prior  opinion).  For  example,  the  optical  contrast  between  tumor 
and  normal  breast  tissue  is  around  3-451  which  gives  the  expected  standard  deviation 
(cqt_Mo)  in  the  optical  properties,  which  can  be  used  to  compute  variance.  The  calibration 
of  the  experimental  data  is  capable  of  giving  a  very  good  estimate  of  normal  tissue  optical 
properties.45,46  Note  that  weight  matrix  containing  the  expected  variance  will  not  impose 
a  hard  constrain  on  the  expected  optical  properties,  but  discourages  update  values  (Ap) 
which  are  above  these  expected  deviation  in  a  given  iteration. 


The  advantages  of  the  GLS  approach  are  that: 

•  It  accounts  for  covariance  among  the  parameters  and  data  points. 

•  It  also  allows  the  individual  data  points/parameters  to  have  different  noise  character¬ 
istics  (variances) 

•  It  constrains  the  problem  through  the  weight  matrices,  to  produce  stable  solutions. 
The  limitations  of  the  procedure  are: 
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•  It  requires  prior  knowledge  about  the  noise  characteristics  of  parameter  and  data  spaces 

•  The  weight  matrices  may  necessitate  computation  of  the  inverse  of  covariance  matrices 
(increasing  run  time  and  memory  requirements) 

•  It  can  generate  unstable  solutions  when  unreasonable  constraints  are  inadvertently 
applied 


With  Spatial  Priors 

Overall,  the  LS  minimization  schemes  using  spatial  priors  can  be  broadly  classified  into  two 
approaches.  (1)  Soft-Priors  (2)  Hard-Priors.  The  following  two  subsections  will  discuss  these 
two  approaches. 

3.D.  Soft- Priors 

In  this  approach,  the  regularization  matrix  L  in  the  Tikhonov  approach  (Eq.  18)  encodes 
the  spatial  information.21,26  Previous  results  have  shown  that  using  the  spatial  priors  in 
this  fashion  does  not  bias  the  image  estimate  when  the  prior  information  is  imperfect.26 
Typically,  the  conventional  image  is  segmented  into  different  regions  depending  on  tissue- 
type  to  generate  the  spatial  constraints.  The  L-matrix  relates  each  nodal  optical  property  in 
the  numerical  model  to  the  other  nodes  in  that  region.26  Two  possible  forms  are  indicated 
below. 


3.D.I.  Laplacian  form21 

!0  if  i  and  j  are  not  in  the  same  region 

—  1/N  if  i  and  j  are  in  the  same  region 

1  if  i  =  j 

where  N  is  the  number  of  sampling  points  (e.g.  nodes  in  a  FEM  mesh)  in  that  region 


(30) 


3.D.2.  Helmholtz  form26 

{0  if  i  and  j  are  not  in  the  same  region 

~  N+(Kh)2  if  1  an(i  J  are  in  the  same  region  (31) 

1  if  i  =  j 

where  N  is  the  number  of  nodes  in  that  region,  k  —  l/l  with  l  being  the  covariance  length, 
and  h  is  the  distance  between  the  nodes. 
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3.E.  Hard- Priors 

In  the  Hard-Prior  approach,  also  known  as  a  parameter-reduction  technique,  the  number 
of  parameters  to  be  estimated  becomes  the  number  of  regions  segmented  from  the  other 
imaging  modality  (spatial-priors).26  Even  though  the  number  of  parameters  to  be  estimated 
reduces  considerably  (relative  to  soft  priors),  the  problem  can  still  be  ill-posed,2  so  an  LM 
approach  was  used  (Eq.  11)  in  this  case  due  to  its  simplicity.  The  main  advantage  of  the 
method  are: 

•  The  problem  is  over-determined,  which  also  implies  J 1  ,J  is  positive  definite 

•  It  is  computationally  efficient 
The  limitations  include: 

•  The  effect  of  error  or  uncertainty  in  the  spatial  priors  can  be  amplified  by  the  technique. 

•  The  DOT  problem  may  still  be  ill-posed  (and  ill-conditioned)  after  the  constraints  are 
added2 

3.F.  Important  Notes  about  Minimization  schemes 

There  are  additional  important  points  about  these  minimization  schemes. 

•  The  weight  matrices  (W#  and  in  the  GLS  scheme  are  computed  before  the 

iterative  reconstruction  procedure  begins  and  are  invariant  during  the  iterative  process. 
The  same  is  true  of  the  Soft-Priors  L-matrix  calculations. 

•  The  first-order  conditions  (Eqs.  5,  15,  and  21)  derived  by  minimizing  the  objective 
functions  (Eqs.  4,  12,  and  20)  in  all  minimization  schemes  appear  on  the  right  hand 
side  (R.H.S)  of  the  update  equations  (Eqs.  11,  18,  and  22)  which  means  that  only  when 
the  R.H.S.  has  reached  zero,  the  solution  reached  the  global  minimum. 

•  Computation  of  weight  matrices,  L-matrices  and  the  Tikhonov  regularization  param¬ 
eter,  requires  a  prior  opinion  about  the  variances  of  the  parameters  and  data.  Here, 
only  the  best  prior  estimates  (BPE)  are  used,  which  means  that  the  actual  variances 
of  the  parameter  and  data  spaces  are  used  in  the  reconstruction  procedure.  Variation 
from  the  best  prior  values  can  be  examined  also,  to  observe  the  effect  of  priors,  but 
that  work  is  beyond  the  scope  of  the  present  paper. 

•  When  spatial-priors  are  used  in  this  study  (as  well  as  in  most  studies),  it  is  assumed 
that  they  are  perfect.  The  effect  of  spatial  prior  uncertainty  on  DOT  inverse  problem 
is  discussed  in  Ref., 23,24,26  and  is  the  subject  of  ongoing  study. 


13 


•  The  covariance  lengths  associated  in  the  weight  matrix  (GLS-AC  form,  Eq.  28),  and 
the  L-matrix  (Helmholtz  form,  Eq.  31)  calculations  are  chosen  to  be  10mm  and  5mm 
respectively.  The  effect  of  covariance  length  on  the  image  reconstruction  is  discussed 
in  Ref.26 

•  In  the  LM  approach  (Eq.  11),  the  Jacobian  is  normalized  by  its  optical  properties.  Also 
a  was  chosen  initially  to  be  1  and  it  was  reduced  by  a  factor  of  10°'25  at  every  iteration 
and  multiplied  by  the  maximum  of  the  diagonal  values  of  JTJ.  The  normalization 
procedure  is  described  in  Ref.52  Moreover,  eight  iterations  were  chosen  for  all  the  LM 
reconstructions,  as  it  has  been  shown  in  the  literature  that  after  this  iteration,  error 
in  the  optical  properties  increases  for  this  particular  problem  and  algorithm.53,54  This 
inherent  instability  can  be  attributed  to  the  fact  that  J 1  J  is  not  positive  definite  in 
DOT. 

•  For  simplicity,  all  the  reconstruction  algorithms  are  tested  only  in  the  two-dimensional 
case.  Comparison  of  three-dimensional  reconstructions  are  left  for  future  investigations. 

3.F.I.  Special  cases  of  GLS  minimization 

The  update  equation  for  the  GLS  scheme,  Eq.  22,  turns  into  the  Tikhonov  case  (Eq.  18)  when 
Wg  =  I  and  WIL_n0  =  A LTL.  Moreover,  if  one  assumes  that  A/i  =  /i  —  /i0,  which  is  equivalent 
to  taking  a  single  step  in  the  iterative  procedure,  then  Eq.  19  maps  into  Eq.  11  with  a  =  2 A. 
Hence,  the  LM  technique  can  be  viewed  as  a  special  case  of  the  Tikhonov  method,  which 
itself  is  a  special  case  of  the  GLS  approach.  It  is  important  however  to  differentiate  LM  from 
the  single-step  Tikhonov  approach  because  LM  requires  a  to  reach  zero  asymptotically  with 
number  of  iterations,  whereas  in  the  Tikhonov  scheme,  A  is  constant.  Moreover,  LM  does 
not  involve  parameters  in  the  objective  function. 

3.F.2.  Stopping  Criterion 

The  importance  of  the  stopping  criterion  in  an  iterative  procedure  can  not  be  ignored.  The 
stopping  criterion  used  in  this  work  is  based  on  the  first-order  conditions  and  data-modcl 
misfit,  which  in  the  limit  ensures  that  the  problem  has  reached  the  global  minima.  The 
iterative  procedure  is  stopped  when  the  L2-norm  of  the  data-model  misfit  (J)  does  not 
improve  by  more  than  10-1°%  or  the  L2-Norm  of  the  first  order  conditions  is  less  than 
1(T17%.  Beyond  these  values,  the  round-off  error  dominates.  This  stopping  criterion  is  more 
robust  because  it  involves  first-order  conditions  as  well. 
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4.  Test  Problem 


This  section  provides  the  details  of  the  test  problem  considered  here.  The  optical  property 
distributions  used  for  the  synthetic  data  (y,  noise  added)  generation  are  shown  in  Fig.  2. 
The  diameter  of  the  domain  was  86  mm.  The  background  optical  properties  were  ya  = 
0.01  mm-1  and  y's  =  1.0  mm-1.  There  were  two  irregular  shaped  targets,  one  in  ya  with  a 
contrast  of  2:1  to  background  and  one  in  y's  with  a  contrast  of  3:1  relative  to  the  background. 
A  mesh  consisting  of  4617  nodes  (corresponding  to  9040  linear  triangular  elements)  was 
used  for  the  generation  of  data.  Sixteen  light  collection/dclivery  fibers  were  arranged  equally 
spaced  on  the  boundary  of  the  circle,  where  one  fiber  was  used  as  the  source  while  all  other 
fibers  served  as  detectors  in  turn  which  produced  a  total  of  240  measurements  (that  is  240 
ln(A)  data  points  and  240  9  data  points).  The  source  was  modeled  as  a  Gaussian  profile 
with  a  Full  Width  Half  Maximum  (FWHM)  of  3  mm  to  represent  the  light  applied,55  and 
was  placed  at  a  depth  of  one  transport  scattering  distance  from  the  tissue  boundary.56 
Noise  levels  of  1%,  3%,  5%  and  10%  were  added  to  the  modeled  data  ([ln(A);  9 ])  to  form 
the  experimental  data  ( y ).  At  the  same  time,  the  variances  in  the  data  were  also  computed 
to  be  used  in  the  reconstruction  algorithms. 

The  actual  reconstructions  and  forward  modeled  data  computation  were  performed  on 
different  FEM  meshes.5'  This  mesh  has  the  same  diameter  (86  mm)  with  1785  FEM  nodes, 
which  corresponded  to  3418  linear  triangular  elements.56  The  expected  distribution  of  optical 
properties  is  given  in  Fig.  3(a)  (first  column).  Background  optical  properties  were  used  as 
initial  estimates  (yo)  hr  the  evaluation  of  reconstruction  methods.  The  number  of  parameters 
to  be  estimated  was  3570  (1785  in  ya  and  1785  in  y's).  The  number  of  data  points  available 
for  reconstruction  was  480  (240  of  ln(A)  and  240  of  9).  The  dimension  of  J  was  480x3570,  W$ 
was  480x480,  and  H//i-/i0  was  3570x3570.  Optical  property  distributions  were  reconstructed 
from  the  data  without  noise  (bias  calculations)  as  well  as  with  noise  levels  of  1%,  3%,  5%  and 
10%.  The  reconstructions  are  repeated  for  the  case  of  3%  noise  in  the  data  with  increasing 
complexity  (targets)in  the  optical  property  distributions. 

5.  Results  and  discussion 

Initially  all  reconstruction  techniques  were  executed  on  a  data-set  without  noise  to  estimate 
the  bias.  Note  that  for  these  calculations  the  variance  was  found  between  the  data  generated 
using  meshes  (described  in  Section-4,  Fig.  2)  with  4617  nodes  and  1785  nodes.  The  results 
without  employing  spatial  prior  information  from  the  reconstruction  techniques  are  given 
in  Fig.  3(a).  The  first  column  shows  the  expected  distribution  for  the  1785  node  mesh 
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used  in  the  reconstruction  and  forward  model  calculations.  The  Tikhonov  approach  failed 
to  recover  the  contrast.  This  was  primarily  due  to  the  choice  of  A,  which  was  based  on 
the  maximum  variance  value,  which  biases  the  problem  to  data  points  that  are  above  the 
average  noise  level.  Since  DOT  is  known  to  have  a  large  dynamic  range  in  the  data  (at  least 
8  orders  of  magnitude55),  this  choice  of  A  deemphasize  the  data  points  that  have  low  or 
intermediate  variance  values.  The  Root-Mean-Squared  (RMS)  errors  between  the  expected 
and  reconstructed  optical  properties  are  plotted  in  Fig.  6.  The  mean  and  standard  deviation 
in  the  reconstructed  images  for  different  regions  (labeled  in  first  column  of  Fig.  3(a))  using 
the  reconstruction  techniques  discussed  till  now  are  given  in  Table-1.  In  the  case  of  no 
spatial  priors,  Levenberg-Marquardt  (LM)  gives  less  bias  in  /ia,  where  as  GLS  performs 
better  in  /% .  The  bias  calculations  were  repeated  with  spatial-priors  and  the  reconstruction 
results  are  presented  in  Fig.  3(b).  These  RMS  errors  in  the  optical  properties  are  also 
plotted  in  Fig.  6.  Surprisingly  the  Soft-Prior  approach  (Laplacian  and  Helmholtz)  performed 
better  than  the  Hard-Prior  strategy.  It  can  also  be  observed  from  Fig.  6  and  Table-1  that 
the  usage  of  spatial-priors  reduces  the  bias  by  at  least  a  factor  of  2. 

Figure  4(a)  shows  reconstruction  results  using  data  with  5%  noise  in  amplitude  without 
employing  spatial  priors.  Once  again  the  Tikhonov  approach  fails  to  recover  the  contrast. 
The  LM  results  are  dominated  by  boundary  artifacts.  Fig.  4(b)  presents  the  results  from 
the  same  data  set  when  spatial  priors  were  employed.  Fig.  5(a)  and  5(b)  show  similar 
kinds  of  effort  for  the  case  of  data  with  10%  noise.  The  RMS  error  in  the  reconstructed 
Ha  and  /i's  images  are  plotted  in  Fig.  6  as  a  function  of  increasing  noise  level.  The  RMS 
error  using  the  LM  approach  increases  with  increasing  noise.  GLS  techniques  perform  very 
well  even  in  the  case  of  10%  noise  (Fig.  5(a)  and  6).  Among  the  GLS  methods,  usage  of 
an  analytical  covariance  form  gives  better  results  (~  13%  less  RMS  error)  in  /ia  and  the 
local  Laplacian  form  performs  slightly  better  (^  3%  less  RMS  error)  in  /%.  In  the  case  of 
employment  of  spatial-priors,  it  can  clearly  be  seen  (from  Fig.  4(b),  5(b)  &  6  and  Table-1) 
that  Hard-Priors  perform  better  in  /i(  reconstruction  when  the  noise  level  is  below  10%. 
Among  the  soft-prior  results,  for  /ia ,  the  RMS  error  linearly  increases  with  increasing  noise 
level  in  the  Laplacian  case  (Fig.  6).  In  /%  reconstructions,  the  performance  of  Laplacian 
and  Helmholtz  are  comparable,  clearly  Helmholtz  performs  slightly  better  (^  5%)  when  the 
noise  level  is  above  3%.  Interestingly,  the  Helmholtz  regularization  emerges  with  the  lowest 
RMS  error  in  /j,a  reconstruction.  This  is  primarily  because  of  the  covariance  length  factor  in 
the  Helmholtz  form  of  the  regularization  matrix  (Eq.  31),  which  ensures  that  the  optical 
properties  covary  within  that  correlation  length  (in  here  it  is  5  mm).  The  same  explanation 
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is  true  for  the  GLS- Analytical  covariance  form  (Eq.  28),  which  performs  better  in  /ia 
estimation.  It  is  also  important  to  note  that  in  the  case  of  a  limited  number  of  wavelengths, 
Srinivasan  et  al58  have  shown  that  5%  error  in  the  optical  property  estimate  (/J,a  and  /i()  can 
lead  to  approximately  45%  error  in  spectral  properties  (Hemoglobin,  Water  Concentrations, 
Oxygen  Saturation,  and  Scattering  Estimates)  of  tissue.  Any  small  improvement  in  the 
optical  property  estimates  would  be  important  for  improvement  in  the  utility  of  this  type 
of  imaging  under  practical  conditions. 

To  emphasize  the  effects  of  complexity  on  the  reconstruction  procedures,  a  set  of 
simulations  were  performed  with  an  increasing  number  of  targets.  Each  target  was  chosen 
to  be  circular  with  a  diameter  of  10  mm.  The  contrast  to  background  optical  properties  was 
2:1.  The  target  locations  and  corresponding  optical  properties  are  shown  in  the  first  column 
of  Fig.  7(a).  The  targets  were  also  labeled  from  1  to  4  (background  is  labeled  as  0).  The 
data  used  in  this  case  had  a  noise  level  of  3%.  A  total  of  4  different  reconstructions  were 
performed  by  adding  each  target  at  a  time  (from  1  to  4).  The  result  of  the  4  target  case  is 
shown  in  Fig.  7.  Corresponding  mean  and  standard  deviation  of  the  reconstructed  optical 
properties  for  different  regions  (labeled  in  Erst  column  of  Fig.  7(a))  are  given  in  Table-2. 
Fig.  8  contains  a  plot  of  RMS  error  in  the  reconstructed  optical  properties  with  increasing 
number  of  targets.  The  RMS  error  increases  with  increasing  number  of  targets  for  every 
reconstruction  algorithm.  Note  that  targets  3  and  4  were  placed  close  to  the  center  of  the 
domain,  where  the  sensitivity  is  low  compared  to  the  periphery.56  Moreover,  increasing 
the  /ra  targets  (from  1  to  2,  target  numbers  1  and  3),  caused  the  RMS  error  to  increase 
by  at  least  30%.  The  same  is  true  with  the  fi’s  targets.  In  the  case  of  multiple  targets, 
the  Helmholtz- type  of  regularization  matrix  resulted  in  the  least  error  in  both  /ia  and  fi's. 
Even  though  the  Hard-Prior  case  performs  very  well  in  terms  of  lowest  RMS  error  for  a 
single  target,  as  the  complexity  (or  number  of  parameters  to  be  estimated)  of  the  problem 
increases,  it  clearly  performs  poorer  than  most  of  the  techniques  presented  (Fig.  8). 

Even  though  the  choice  of  Tikhonov  regularization  parameter  (A)  given  by  Eq.  14  is 
the  optimal,  the  other  common  way  is  to  use  L-curve.59  The  L-curve  for  DOT  is  much 
shallower60  similar  to  the  estimation  problem  in  Electrical  Impedance  Tomography  (EIT), 
which  poses  a  problem  in  selection  of  A  and  is  shown  being  unreliable  in  Ref.5' 

Table-3  gives  the  computational  time  per  iteration  for  each  of  the  reconstruction  technique 
(in  these  two-dimensional  cases)  on  Pentium  IV  (dual  core)  2.8  GHz,  2GB  RAM  Linux  work 
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station.  GLS  schemes  take  little  more  computation  time  than  the  Tikhonov  minimization, 
as  expected  Hard-Priors  took  the  least  computation  time. 

Overall,  the  inclusion  of  spatial-priors  has  an  important  positive  effect.  The  errors  in  the 
estimated  optical  properties  are  also  reduced  by  at  least  a  factor  of  2  with  spatial-information. 
The  reconstructed  images  also  contain  the  fine  features  extracted  from  conventional  imaging 
modalities.  Through  the  incorporation  of  the  individual  variability  of  the  data  points  and 
optical  parameters  (GLS  scheme),  reconstruction  performs  better  even  when  the  noise  level 
in  the  data  is  high.  It  is  also  important  to  note  that,  as  mentioned  before,  iteration  number  8 
(which  is  the  best  result  in  terms  of  lowest  RMS  error)  is  chosen  for  RMS  error  calculations 
in  LM  approach,  after  this  iteration,  the  solution  becomes  unstable.  Whereas  the  rest  of  the 
approaches  yield  stable  solutions  (error  in  optical  properties  did  not  increase  with  increas¬ 
ing  iterations).  When  the  individual  data  point  variances  were  not  considered  (Tikhonov 
approach),  the  reconstruction  algorithm  may  not  have  the  ability  to  recover  the  contrast  in 
the  target.  Moreover,  simultaneous  estimation  of  both  absorption  and  scattering  coefficients 
causes  cross-talk  between  the  two  parameter  estimates.  Even  with  error-free  spatial-priors, 
as  the  complexity  of  the  estimation  problem  (or  number  of  targets)  increased  for  a  given 
noise  level  in  the  data,  the  parameter-reduction  (Hard-Priors)  technique  failed  to  give  the 
best  estimates  of  the  optical  properties  clue  to  its  LM  implementation. 

6.  Conclusions 

The  diffuse  optical  tomography  inverse  problem  is  often  solved  by  Levenberg- 
Marquardt /modified  Tikhonov  minimization.  A  generalized  approach  for  diffuse  optical 
tomographic  imaging  which  incorporates  the  expected  variability  of  the  data  noise  and 
magnitude  of  the  optical  parameter  variation  is  presented  as  a  structured  weight-matrix 
regularization.  It  is  also  shown  that  Tikhonov  minimization  and  the  Levenberg-Marquardt 
approach  are  special  cases  of  this  generalized  Least-Squares  (GLS)  minimization  formalism. 
Weight-matrices  that  are  used  in  this  reconstruction  procedure,  consisting  of  the  variance 
and  covariance  among  the  data  points  and  optical  properties,  penalize  the  solution  to  match 
the  modeled  data  with  the  experimental  data  more  appropriately.  This  frame-work  can 
also  be  used  to  incorporate  structural  information,  given  by  MR,  CT  or  other  imaging 
modalities  when  the  two  are  acquired  on  the  same  tissue  volume.  Using  a  test  problem, 
all  of  these  techniques  are  studied  in  terms  of  data  noise  level  and  test  field  complexity 
and  a  uniform  comparison  was  made  using  the  same  implementation  scheme  for  each 
minimization  method.  Even  with  highly  noisy  data,  the  GLS  approach  gives  meaningful 
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reconstruction  results.  It  appears  that  the  standard  Levenberg-Marquardt  approach  may 
be  unstable  for  the  DOT  problem.  It  is  also  shown  that  consideration  of  the  individual 
variances  of  data-points  is  the  key  for  an  estimation  procedure  to  recover  high  optical 
contrast.  Employing  spatial  information  reduced  the  errors  in  the  reconstruction  results 
by  at  least  a  factor  of  2.  Parameter-reduction  using  spatial-priors  can  produce  erroneous 
results  when  the  noise  level  is  high.  The  same  is  true  for  increasing  numbers  of  targets. 
Future  work  includes  investigating  various  approaches  for  incorporating  spatial-priors 
into  the  GLS  scheme  with  experimental  data  sets.  Moreover,  a  thorough  examination  of 
these  techniques  in  three-dimensional  case  will  be  taken  up  as  a  future  investigation.  The 
computer  algorithms  and  test  data  used  in  this  paper  (along  with  some  additional  material) 
are  given  at  this  web  page.61 


Appendix 

A.l  Terminology 

DOT-Diffuse  Optical  Tomography. 

/ia(r)-Optical  absorption  coefficient  of  tissue. 

/./((r) -Reduced  (or  transport)  scattering  coefficient  of  tissue. 

_D(r)-Optical  diffusion  coefficient  of  tissue  =  y,,  (r)+M'  (r)]  • 

/^-Parameters  (generalized)  to  estimate  =  [Z)(r);/ta(r)]. 

/loTPrior  value  of  the  parameters  (initial  guess,  generally  obtained  from  prior  calibra¬ 
tion  of  data45,46). 

F(/i)-Forward  Model. 

G(/i)-Modelcd  data  (G  -  Sampled  Forward  model  =  S'{F}). 

A-  Amplitude  of  the  signal. 

0-Phase  of  the  signal. 


19 


^/-Measured  data  =  [ln(A);  6  . 


||X||-L2-norm  of  vector  X  =  \/Y2iLi  Xi2 ■ 

5-Data-Model  misfit  =  y  —  G(fi). 

W^-Weight  matrix  for  5  =  ( cov(5 ))”1  (Appendix  A-4  of  Ref.47). 

W/7,t_,t0-Weight  matrix  for  /n-fi0  =  (con(/i  —  /io))^1  (Appendix  A-4  of  Ref.4'). 

A-Tikhonov  regularization  parameter. 

L-Tikhonov  regularization  matrix. 

/-Identity  matrix. 

<r2-Variance 

J-Jacobian  of  the  sampled  forward  model  =  9Gq^  ■ 
fl-Objective  function. 

Error-True  value  -  Estimated  value  (prediction). 

Bias-Difference  between  the  true  optical  property  distribution  and  estimated  optical 
properties  in  the  case  of  model  generated  data  (without  adding  the  noise). 

Ill-posed-Small  changes  in  the  data  can  cause  large  changes  in  the  parameters. 

Ill-conditioned-The  condition  number  (ratio  of  largest  singular  value  to  smallest 
singular  value)  is  large,  which  implies  the  inverse  solution  would  not  be  unique. 

Ill-determined-(or  under- determined)  The  number  of  independent  equations  are 
smaller  than  number  of  unknowns. 
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Unstability-Error  gets  amplified  with  iterations. 


LM-Levenberg-Marquardt  minimization  (Sec.  3. A). 

Tikhonov-Tikhonov  minimization  scheme  without  spatial-priors,  L  =  I  (Sec.  3.B). 

GLS- AC-Generalized  least  squares  minimization  scheme  (Sec.  3.C)  with  analytical 
covariance  form  for  (Eq.  28). 

GLS-LL-Generalized  least  squares  minimization  scheme  (Sec.  3.C)  with  local  Lapla- 
cian  form  for  (Eq.  29). 

Laplacian-Tikhonov  minimization  scheme  in  the  case  of  soft-priors  (Sec.  3.D)  where 
L  approximates  Laplacian  form,  defined  by  Eq.  30. 

Helmholtz-Tikhonov  minimization  scheme  in  the  case  of  soft-priors  (Sec.  3.D)  where  L 
approximates  Helmholtz  form,  defined  by  Eq.  31. 

Hard-Priors-Parameter- reduction  technique  based  on  spatial  priors  (Sec.  3.E). 


A. 2  Tikhonov  Regularization  -  Singular  Values 

It  is  interesting  to  examine  Tikhonov  regularization  from  the  point  of  view  of  singular 
values.  If  one  rewrites  the  update  equation  (Eq.  19)  as 


J  +  AI 


A/q 


JT*-1  +  c 


(32) 


where  C  =  A(/q_i  —  go),  as  it  is  a  constant  vector  for  a  chosen  iteration  i.  Singular- Value 
decomposition  (SVD)  of  J  gives 

J  =  USVT  (33) 

where  U  and  V  are  orthonormal  matrices  containing  the  singular  vectors  of  J,  i.e.  UTU  =  I 
and  VTV  =  I.  S  is  a  diagonal  matrix  containing  the  singular  values  (Si)  of  J.  Substituting 
this  into  update  equation  (Eq.  32)  generates 


VSTUTUSVr  +  XI  A =  V&WSi-!  +  C 


TttT> 


(34) 
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Using  the  orthonormal  properties  of  U  and  left  multiplying  by  VT  on  both  sides  of  Eq.  34 
yields 


V 1  VS 1  SV 1  +  XV1 


A  m  =  V1  VSTUT8i-i  +  VTC 


(35) 


Now  using  the  orthonormal  properties  of  V  and  rearranging  the  terms  leads  to 


STS  +  XI  VTA ^  =  STUTdi_ i  +  VTC 


(36) 


Taking  the  inverse,  left  multiplying  by  V  and  simplifying  the  result  gives 


A  m 


V 


i-i 


STS  +  XI  '  STUT5i_1  +  V'1'C 


(37) 


Writing  Eq.  37  in  the  form 


Aim  =  VDP 


(38) 


where  P 
form 


STUT5i-\  +  VTC  ,  a  column  vector,  and  D  is  a  diagonal  matrix  which  has  the 


D(i,j) 


0  if  i^j 
A~x  'fi=3 


(39) 


Similar  expressions  hold  for  L  J63  in  Eq.  18.  Considering  the  case  A  =  0,  one  can  clearly 
see  that  for  an  ill-conditioned  matrix  J,  implying  some  of  the  singular  values  are  almost 
zero  (Si  ~  0),  the  inversion  becomes  unstable  (some  of  the  diagonal  values  of  D  become 
infinite).  By  using  Tikhonov  regularization,  even  when  Si  =  0,  the  inversion  procedure  is 
stabilized  (Eq.  39).  The  A  act  as  a  filtering  factor,  giving  the  name  Tikhonov  filtering 63  for 
this  procedure.  Moreover,  as  this  A  damps  the  amplification  of  the  diagonal  values  of  D  for 
smaller  values  of  St  in  Eq.  39  ,  this  is  also  known  as  damped  least  squares  minimization 
procedure.63 
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Figure  4:  Reconstruction  results  (top  of  the  first  row,  abbreviations  are  given  in  the 
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spatial  priors.  The  top  row  gives  images  of  pa  and  bottom  row  shows  fi's  images. 

Figure  6:  A  plot  of  the  RMS  error  in  the  estimated  optical  properties  is  shown  as  a 
function  of  increasing  noise  level  for  all  reconstruction  techniques. 

Figure  7:  Reconstruction  results  (top  of  the  first  row,  abbreviations  are  given  in  the 
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Figure  8:  Plot  of  the  RMS  error  in  the  estimated  optical  properties  is  shown  for  in¬ 
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Tables 


Methods 

Noise  level 

Region- 0 

Region- 1 

Region-2 

Actual 

- 

0.01 

0.02 

0.01 

0% 

0.0101±0.001 

0.0172A0.0023 

0.0105A0.0005 

LM 

5% 

0.0102A0.0016 

0.0125A0.0016 

0.0123A0.0011 

10% 

0.0103A0.0029 

0.0132A0.0026 

0.0118A0.0023 

0% 

0.0102A0.0005 

0.0117±0.0003 

0.0117±0.0002 

Tikhonov 

5% 

0.0102A0.0004 

0.0114A0.0002 

0.0112A0.0001 

10% 

0.0102A0.0003 

0.0108A0.0009 

0.0107A0.0005 

0% 

0.01±0.001 

0.015±0.0011 

0.0112A0.0003 

GLS-AC 

5% 

0.0101A0.0014 

0.0146A0.0012 

0.0106A0.0004 

10% 

0.0101A0.0013 

0.0136A0.0009 

0.0111±0.0008 

0% 

0.01±0.001 

0.0152A0.0012 

0.0113A0.0003 

GLS-LL 

5% 

0.0101A0.0016 

0.0149A0.0015 

0.0108A0.0006 

10% 

0.0101A0.0016 

0.0138A0.0009 

0.0112A0.001 

0% 

0.0098A0.0001 

0.0212A0.0001 

0.0112±0.0001 

Laplacian 

5% 

0.0098A0.0002 

0.0247A0.0001 

0.0097A0.0001 

10% 

0.0095A0.0001 

0.0276A0.0002 

0.0157A0.0128 

0% 

0.0099A0.0001 

0.019±0.0002 

O.OllliO.OOOl 

Helmholtz 

5% 

0.0099A0.0002 

0.0193A0.0002 

0.0099A0.0001 

10% 

0.0098A0.0002 

0.0174±0.0002 

0.0136A0.0001 

0% 

0.0099 

0.0218 

0.0116 

Hard-Priors 

5% 

0.0098 

0.0218 

0.0131 

10% 

0.0098 

0.018 

0.0166 

Table  1(a) 
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Methods 

Noise  level 

Region- 0 

Region- 1 

Region-2 

Actual 

- 

1.0 

1.0 

3.0 

0% 

1.0356T0.2364 

0.9995±0.0359 

2.3758±0.5160 

LM 

5% 

1.075±0.0357 

1.0555±0.3254 

1.8215±0.3144 

10% 

1.2672±0.9086 

1.3111±0.4128 

1.7111±0.6112 

0% 

1.0096±0.0397 

1.1153±0.0260 

1.1644±0.0251 

Tikhonov 

5% 

1.0111±0.0004 

1.0912±0.0189 

1.0934±0.0104 

10% 

1.0107±0.0216 

1.0441±0.0062 

1.0416±0.0035 

0% 

1.0034±0.0688 

1.0335±0.0199 

1.6838±0.1961 

GLS-AC 

5% 

1.0008A0.0916 

1.0670±0.0362 

1.6972±0.2037 

10% 

0.9987A0.0831 

1.0761±0.0343 

1.3703±0.0773 

0% 

1.0022A0.0693 

1.03±0.0183 

1.7801±0.2573 

GLS-LL 

5% 

0.9998A0.1035 

1.0567±0.0329 

1.8502±0.3034 

10% 

0.9981±0.0947 

1.0839±0.0425 

1.4271±0.0990 

0% 

0.9918±0.0155 

0.9429±0.0015 

2.8207±0.0491 

Laplacian 

5% 

0.9895±0.0202 

0.8559±0.0036 

3.6931±0.1551 

10% 

1.0103±0.0124 

0.7447±0.0011 

1.9884±0.0096 

0% 

0.9878±0.0154 

1.0518±0.0018 

2.7833±0.0854 

Helmholtz 

5% 

0.9813A0.0199 

1.1204±0.0081 

3.4252A0.1947 

10% 

0.9884A0.0121 

1.2766A0.01 

2.1761±0.0382 

0% 

0.9919 

0.9266 

2.7332 

Hard-Priors 

5% 

0.9874 

1.0358 

2.345 

10% 

0.9854 

1.3899 

1.822 

Table  1(b) 
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Methods 

Region- 0 

Region- 1 

Region-2 

Region-3 

Region ~4 

Actual 

0.01 

0.02 

0.01 

0.02 

0.01 

LM 

0.0101±0.0004 

0.0113±0.0001 

0.0112±0.0002 

0.0111±0.0003 

0.011±0.0002 

Tikhonov 

0.0102±0.0004 

0.011±0.0001 

0.0112±0.0001 

0.0109±0.0001 

0.011±0.0001 

GLS-AC 

0.0102±0.0009 

0.0129±0.0003 

0.0111±0.0003 

0.0114±0.0003 

0.0113±0.0003 

GLS-LL 

0.0102±0.0011 

0.0133±0.0004 

0.0115±0.0004 

0.0113±0.0003 

0.0113±0.0002 

Laplacian 

0.01±0.0002 

0.0181±0.0001 

0.0105±0.0001 

0.0152±0.0001 

0.0158±0.0001 

Helmholtz 

0.01±0.0002 

0.0169±0.0001 

0.0115±0.0001 

0.0149±0.0001 

0.0158±0.0001 

Hard-Priors 

0.01 

0.0158 

0.0126 

0.014 

0.0158 

(a) 


Methods 

Region- 0 

Region- 1 

Region-2 

Region-3 

Region-^ 

Actual 

1.0 

1.0 

2.0 

1.0 

2.0 

LM 

1.0063±0.0986 

1.1333±0.0027 

1.24±0.0623 

1.1191±0.0396 

1.097±0.0366 

Tikhonov 

1.0051±0.0217 

1.0341±0.0019 

1.0575±0.0073 

1.0321±0.0056 

1.0329±0.0043 

GLS-AC 

0.9993±0.0489 

0.9885±0.0139 

1.2486±0.0447 

1.021±0.0234 

1.1184±0.0076 

GLS-LL 

0.9987±0.0553 

0.9764±0.0127 

1.2726±0.0596 

1.0271±0.0262 

1.1422±0.0105 

Laplacian 

0.9886±0.0163 

1.0891±0.0023 

2.3799±0.0242 

1.3445±0.0043 

1.4044±0.0036 

Helmholtz 

0.9899±0.0164 

1.1499±0.0037 

2.1122±0.0386 

1.3382±0.0079 

1.3521±0.0066 

Hard-Priors 

0.9856 

1.3712 

1.7319 

1.4471 

1.5255 

(b) 


Table  2 
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Reconstruction  Method 

Computation  time  per  iteration 

LM 

17.92  Sec 

Tikhonov 

21.28  Sec 

GLS-AC 

23.39  Sec 

GLS-LL 

23.39  Sec 

Laplacian 

22.78  Sec 

Helmholtz 

22.78  Sec 

Hard-Priors 

10.73  Sec 

Table  3 
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Figure  1 
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Figure  3(a) 
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Figure  3(b) 
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Figure  4(a) 
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Figure  4(b) 
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Figure  5(a) 
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Figure  5(b) 
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Figure  7(a) 
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Figure  7(b) 
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Figure  8 
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